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ABSTRACT 


With recent delivery of littoral combat ships, the impact of operating in shallow or littoral 
ocean domain (LOD) during the duration of their life cycle is of interest, and a shock trial 
or hardening test and validation for this class is needed. For this study, the theories of 
underwater shock phenomena as applied within the boundaries of LOD, specific to the 
Eulerian fluid domain were conducted. The results of varying ocean depth show clear 
distinction in UNDEX characterization at depths shallower than 300ft. Varying charge 
size and depth showed that charge size of less than 3001bs of HBX-1 displayed a linear 
relationship while changing the charge depths to near water-air or water-bottom interface 
also resulted in amplified characteristics of UNDEX parameters. In addition, varying 
lateral boundary showed that as its distance is brought inside the radius of bulk 
cavitation, the UNDEX behavior also became increasingly chaotic due to similar effects 
seen in the shallower bottom depth. Lastly, adding blocked cells prior to a full scale 
coupled run showed that fluid behaves more erratically as these small rigid boundaries 
are situated within the radius of forming bulk cavitation. 
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I. 


INTRODUCTION 


A. BACKGROUND 

The theory of underwater explosion (UNDEX) phenomena has been studied 
extensively since post World War II. For United States Navy (USN), its motivation was 
to analyze damages incurred on surface and subsurface vessels during direct impact and 
near misses from underwater weapons. Since then, advancement in computing power, 
extensive understanding of the UNDEX theory and theory of finite element have led to 
solidifying efforts in validating shock trial analysis of newly commissioned ships via 
modeling and simulation prior to an actual at sea shock trial analysis. The result is 
software that utilizes UNDEX theory, based on theoretical works of Arons, Cole, Snay 
and Taylor, called Dynamic System Mechanics Advanced Simulation (DYSMAS) [1]. 
DYSMAS is a result of the combined efforts of Gennan Defense Contractor, 
Industrieanlagen-Betriebsgesellschaft mbH (IABG) and Naval Surface Warfare Center 
Indian Head Division (NSWC - IHD), as well as, Lawrence-Livermore, led by the USN. 
It is a six degree of freedom (DOF) finite element hydrocode used to model and simulate 
fully coupled, fluid-structure interaction (FSA) problems of an UNDEX event on a 
surface or subsurface vessel. 

With the recent delivery of the littoral combat ship (LCS), the impacts of 
operating in a shallow or littoral ocean domain (LOD) during the duration of the ship’s 
life cycle are of interest. Thus, a shock trial or hardening test and validation for this class 
is needed. However, due to current Department of Defense’s (DOD) budgetary 
constraints and rising LCS class issues, an actual shock trial analysis for LCS class has 
been postponed indefinitely [2], Nevertheless, continued study of FSA within the LOD 
during UNDEX is needed in order to verify the sea and battle worthiness of LCS class. 

From Naval Postgraduate School’s (NPS) Shock and Vibration Computation Lab 
(SVCL), several studies have already been conducted to address the factors of UNDEX 
simulations of littoral or shallow water by Walters, Arbogast, and Santiago regarding 
explicitly modeled solid Lagrangian ocean floor, influence of shock-induced air bubble 
collapse and cavitation effect on surface ship model, respectively. As Walters suggests, 
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the distinct advantages of using Lagrangian solid bottom was the ability to model 
contoured bottom shapes and for the various solid bottom types created, there were 
significant differences between the initial bottom reflections for the different contours 
[1]. He further suggests that the most important bottom contour effect was the distortion 
to the gas bubble and its associated first pulse timing. Santiago’s study of bulk cavitation 
(BC) effect on fluid-structure interaction (FSA) shows that the effects of cavitation 
closure can be quite significant and inflict another impulsive force upon the ship as 
reported empirically during an actual shock trial report of DDG 53, in addition to the 
incident wave, bubble pulsations, bottom reflection, or bottom refractions [3], Lastly, 
Arbogast’s work showed that the bubbles in the vicinity of marine vessels, depending on 
the size and location and proximity of these bubbles, created a buffering effect from 
shock induced by an UNDEX [4]. 

B. SCOPE OF RESEARCH 

The importance of FSA and depth of such UNDEX theory in evaluating and 
analyzing shock design of a naval vessel cannot be understated. However, in order to 
understand the FSA, first, a solid understanding of how fluid behaves during UNDEX 
must be clearly sought out. As Figure 1 shows, there are numerous theoretical and 
experimental analysis that deals with the interdisciplinary studies of Naval Ship Shock 
Design and Analysis [5]. In short, it is primarily divided into structural analysis that 
stems from a thorough understanding of dynamic fluid characteristics and fluid-structure 
interactions during UNDEX through the use of finite element modeling and simulations, 
as well as, real world experimental validations. For this study, the theories of underwater 
shock phenomena as applied within the boundary of shallow water or LOD will be of 
focus. 

Building on the knowledge of previous researches conducted here at NPS and 
various UNDEX studies published over the years, a more in-depth analysis and 
characterization of LOD is attempted. Furthermore, a clear definition of LOD, in terms 
of UNDEX parameters, especially the influence on BC is sought out. In addition, the 
implementation of various boundary conditions for both bottom and side conditions of 
fluid domains were compared and contrasted in the current Eulerian fluid model, gaining 
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better understanding of their effect on shallow water UNDEX phenomena. In short, 
several variations of typical benchmark LOD problems were modeled and simulated 
using DYSMAS to gain better insights to its resulting UNDEX parameters including the 
impacts that BC have on the overall response to fluid behaviors. 
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II. UNDERWATER EXPLOSION THEORY 


A. INTRODUCTION TO CONVENTIONAL, NON-NUCLEAR EXPLOSION 

A conventional, non-nuclear explosion is a chemical reaction in a material or 
substance that converts the original material into high pressure gas and temperature at 
extreme and catastrophic speed [6], In general, any explosive material, regardless of its 
state, is an unstable compound that becomes more stable during the chemical reaction 
that takes place post ignition and detonation. The temperature within the explosion of a 
gas is usually in the order of 3,000° C (5,432° F) and the pressures are greater than 
50,000atm (734,797psi) [6]. A reaction of this magnitude can only be initiated with 
sufficient energy provided to a specific place within the explosive prior to detonation and 
is done usually by means of a fuse. 



Figure 2. Conventional, Non-nuclear Explosion on Land, from [7] 

During detonation process, the developing intense heat and pressure are sufficient 
to set up the explosive material surrounding the fuse. As a result, a very narrow 
boundary between the exploding material and its surrounding develops as a product of 
extremely high temperature, pressure and energy. This distinctly defined rapidly 
advancing discontinuity, kn own as a “detonation” or “blast wave,” travels with a velocity 
of several thousand feet per second and destroys everything in its path (Figure 2) [6]. 
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The way in which blast waves or detonating disturbance propagates depends on the 
physical or chemical properties of explosive materials and the surrounding medium—air, 
water or vacuum—at which the explosion takes place. From the standpoint of 
destructiveness of explosive materials and impact to its surrounding, the process of 
detonation, propagating shock waves (SW) and energy transfer that takes place through 
the means of pressure or velocity changes to its surrounding becomes increasingly 
important. 

B. UNDERWATER EXPLOSION 

Conventional UNDEX is a non-nuclear explosion that takes place at the surface or 
subsurface of water. The characteristics of UNDEX is similar to that of conventional 
explosion discussed previously, except due to water as medium, the initial mass of 
explosive becomes unstable hot mass of gas at tremendous pressure that radiates 
spherically affecting its surrounding (see Figure 3(a)). The energy expenditure during 
UNDEX detonation accounts for predominantly the SW energy (45%) and bubbles (45%) 
with approximately 10% unaccounted for (see Figure 3(b)) [8], Within the energy 
expenditure partition of SW and bubbles, local and BC are also accounted for some of the 
energy that was imparted into its surrounding water during detonation. 

For this research, water is characterized as a homogeneous fluid capable of 
supporting only limited shearing stresses; hence, a medium and its volume of which can 
readjust itself to displacement of boundaries by result of flow [6]. In addition, as 
conservation of energy and work implies the change in its pressure on confined mass 
results in compression or work done to that system, which creates change in volume of 
such mass. Hence, when pressure is applied to a localized region, the compressibility of 
water makes it possible for it to transmit energy as a wave disturbance to other parts 
within the domain with a finite velocity which results in a local motion of water and 
pressure changes. 
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(b) Energy Partition in UNOEX 



Figure 3. (a) Shock Wave and (b) Energy Partition at Time of Detonation, 

from [5], after [8] 

If the pressure applied to localized regions within the water domain is small 
enough, the rate of propagation or wave of disturbances can be independent of its 
pressure magnitude. This propagating SW travels at about 4,900ft/sec (fps), in sea water 
at 18° C (64.4° F) [6]. Of note, this velocity of the propagating disturbance is also 
dependent on other characteristics of the water, like salinity and its density. If this wave 
motion is simplified to one-dimension so that plane waves are assumed, vice spherical, 
the wave travels without significant changes in its characteristics. However, if the 
detonating waves are radiated from a spherical source, there is an inverse relationship 
between the amplitude of the radiating SW and distance from its source. During the 
diverging propagation of spherical wave, motions of its surrounding water are affected 
due to “surge or afterflow” that takes place in reaction to propagating waves [6]. 

As this SW propagates, the effects it has on the surrounding medium and various 
boundaries like water-air and water-bottom interfaces, cause additional tensile or 
compressive wave that induces additional havoc to existing underwater objects. When 
initial SW reaches water-air interface, the differences in tensile and compressive 
strengths of water and air causes massive cavitation called, BC as it spallates the water 

7 








particles on the surface and separates this layer from below. The initial bubble with its 
detonating cycles of its expansion-collapse motion as it rises to the surface and spews out 
exploding material into the atmosphere due to buoyancy, it too creates another 
catastrophic havoc to anything in its path or in the vicinity there of. Hence, the five 
primary events that take place upon detonation during UNDEX are initial SW 
propagation, formation of BC, cyclic bubble phenomena, secondary or third pressure 
pulses due to boundary reflections and surface phenomenon. 

1. The Shock Wave and Acoustic Wave Assumptions 

The initial cause of disturbance to the water in an UNDEX is the arrival of shock 
pressure wave reacting between the boundaries of explosives and water [6]. This 
pressure distribution, usually in the order of 34,0001b/in“ (psi) for a 3001bs of TNT, 
begins to rapidly decay as it propagates followed by intense peak pressure wave in an 
outward motion of the water (figures 4 and 5) as the extremely dense mass of exploding 
gas continues to spherically expand as detonation takes places. As its SW pressure 
diminishes, so too does the pressure in the water fall off rapidly, as shown in three 
different target distances of 5, 50 and 500ft from the point of detonation (Figure 4). 
Lasting for a few milliseconds at most, a discontinuity exists in the initial SW pressure 
rise that is followed by a distinct exponential decay, shown in Figure 5 [6]. 



Figure 4. Pressure Distribution at (a) 5 ft (b) 50ft and (c) 500ft from 

Detonation, from [6] 
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The pressure distributions of a 3001bs TNT charge, shortly after detonation at 
target distance of 50 and 500ft at three instants of time from the center of TNT charge are 
depicted in Figure 5. For comparison, the initial wave in previous Figure 4(a) is also 
depicted in acoustic wave form by dotted lines for (b) 50 and (c) 500ft. The obvious 
scaling law that emerges when different charge sizes are used is called, “principle of 
similarity.” The principle of similarity states that “if the linear size of charge is changed 
by factor of k, the pressure conditions will also change at new distance and time scales k 
times as large as the original ones are used” [6]. 



Figure 5. Shock Wave Pressure Profile at Two Distances from Detonation, 

from [6] 


The radially propagating disturbance as a wave of compression in water called, 
“shock wave” has following characteristics [6]: 

a. The initial velocity of propagation at or near charge and water boundary is 
approximately 5,000fps at detonation and drastically falls to its “acoustic” 
values as wave advances outward (Figure 4) 

b. The pressure level in spherically propagating wave falls off more rapidly 
than with inverse of distance near the charge at detonation, but eventually 
converges to this behavior as it reaches larges distances. 

c. The region of highest peak pressure is at target location closest to the 
charge at detonation and its wave profile broadens gradually as the wave 
propagates outward. Same characteristics can be observed at various target 
location with reduced peak pressure. 
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As mentioned, since the initial SW velocity will be drastically reduced and 
eventually converged to acoustic speed as it propagates, it is assumed that the SWs are 
approximated as acoustic plane waves with small compression and speed. As a result, 
following Equation (2.1) is derived, 



( 2 . 1 ) 


which states that the added fluid pressure (P) from propagating shockwaves is a function 
of water’s mass density (p 0 ), acoustic wave speed (C 0 ) and its propagating velocity (u). 
The medium in which SW travels, water in this case, are considered to be constant 
properties and as mentioned previously, this relationship assumes that water is uniform, 
homogenous, compressible and inviscid. 



Figure 6. UNDEX Geometry and Various Paths of Detonating Shockwave, 

after [5] 

There are basically four primary paths that a propagating shockwave can travel 
during an UNDEX event before reaching its target (Figure 6): (1) Direct Path (2) Surface 
Reflected or Rarefacted (3) Bottom Bounce and (4) Bottom Refracted Waves. 
Depending on the UNDEX geometry of the explosion in terms of charge depth, size and 
ocean depth, these shockwave travels in the order of listed and the nature of applied 
pressure waves are either compressive (+) or tension (-) as experienced by the target due 


10 





to varying densities at water, air or bottom interfaces. Of note, the tensile pressure wave 
of (2) surface cut-off is due to water’s inability withstand tensile force at air-water 
interface; hence, creating a region of expansion in water molecule rather than a 
compressive stress near detonation and other traveling path of pressure wave. 



Figure 7. Arbitrary Interface Geometry, from [1] 

The initial shockwave’s behavior at an interface such as water-air or water-bottom 
(seabed) can be approximated by using combinations of Snell’s law and pressure-velocity 
continuity at the interface [5]. An arbitrary two different medium interface is shown in 
Figure 7. Using Snell’s law, the relationship of various medium’s acoustic speed (C) and 
its angle of incident (a) to the approaching interface can be referenced to its normal 
vector of that interface or summarized in Equation 2.2 below: 


c. 


c. 


c. 


sin a 


sin ol 


( 2 . 2 ) 


Now, if the pressure-velocity continuity is applied, equations (2.1) and (2.2) can 
be combined to describe the interface, yielding the following relationship between the 
pressure and incidental velocity, as well as, reflected, refracted or transmitted 
shockwaves. 
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Therefore, when above relationships are applied to varying boundaries, the 
resulting pressure wave becomes either compressive or tensile waves. At rigid boundary 
(P 2 C 2 » piCi), such as the water-bottom interface, the wave reflections are compressive 
with the magnitude of both incident and reflected waves pressure and velocity being 
equal, since by definition, the velocity of a rigid boundary is zero [1]. For a non-rigid 
boundary condition (P 2 C 2 « P 1 C 1 ), such as water-air interface, the same relationships 
applied as above results in a tensile reflection or rarefaction waves. 


P(t) of Initial Shockwave for W=60lbf of HBX-1 
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Figure 8. MATLAB Figure of Pressure History, from [1] 

The simultaneous passage of a compressive and tensile wave through a point at 
different times causes a unique effect creating a discontinuous drop, is called the pressure 
cutoff shown in pressure-time history in Figure 8. At this point, the abrupt halt to the 
impulse induced on a target at the point of cut-off can also be observed [1]. This 
dependence of UNDEX behaviors on its geometry is depicted in Figure 9 and using such 


Pressure Cut-off 
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geometry, pressure and vertical take-off velocity (VTO) can be calculated to analyze 
target behavior at various instants of time. 



Figure 9. UNDEX Geometry, from [1] 

The obvious interdependence of charge weight and wave propagation distance for 
spherical charge gave rise to empirical formulas to approximate several important 
parameters of UNDEX that have been continuously developed and tested by Cole and 
others since WWII [5]. The two primary formulas are the peak pressure (P max ) and decay 
constants (0) shown in Equations (2.5) and (2.6): 

(2.5) 

( 2 . 6 ) 

where charge weight (W) is in pounds, the wave propagation distance (R) is in feet, the 
peak pressure (P max ) in psi and the decay constant (0) in milliseconds (msec). The 
constants K and A are depended on the explosive materials used and were found 
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experimentally prior to this research through numerous field testing during various 
historical researches [5]. 

2. Typical UNDEX Problem 

A typical UNDEX problem with an explosive charge of weight, W, detonated in 
water at depth, D, is considered in Figure 10 and its corresponding pressure wave profile 
as time passes is shown in Figure 11 [9], The water depth for this problem is assumed to 
be infinitely deep; therefore, negligible when it comes to bottom bounce interaction and 
considerations. In order to measure various parameters of UNDEX event taking place, an 
absolute pressure gage is located at target location of point (X, Y) or 1 to measure the 
pressure changes at this location. As detonation takes place, a compressive SW 
propagates radially or spherically from its source and travels along the line 0 - 1 in 
Figure 10 as it reaches the target or gage location [9]. 



Figure 10. Incident and Reflected Wave Paths, from [9] 

The absolute pressure gage shows constant pressure at its hydrostatic depth 
denoted by line A - B in Figure 11 prior to detonation. This pressure at this depth is sum 
of atmospheric and hydrostatic pressure at the gage depth and upon arrival of SW at 
target point, the pressure jumps to point C, followed by exponential decay along line C - 
D is observed [9]. As will be shown for most of the cases studied in this research, this 
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pressure is at or close to 14.7 psi or latm; for all practical purpose, at or near surface of 
the water. As shown previously in Figure 6 and now Figure 10, the line 0-2 represents 
the compressive path of SW as it travels to the water-air interface. Next, the SW takes 
the path of 2 - 1, where a rarefaction or tensile wave is observed due to difference in 
acoustic impedance of air and water mentioned previously. 



Figure 11. Absolute Pressure at a Point, from [9] 


The arrival of this rarefaction or surface reflected wave at the gage location is 
tenned surface cutoff and is illustrated in Figure 11 by a sudden drop in the absolute 
pressure from point D to three different locations, E, F or G, depending on the vapor 
pressure of the surrounding water [9]. This concept of surface cutoff is made clearer by 
depicting the rarefaction wave as emanating from an image charge, W;, and propagating 
along line 3 - 2 - 1 in Figure 10. This method of including image charge to calculate 
surface cutoff characteristics is known as method of images and is valid since the 
distance 3-2-1 equals the distance 0 - 2 - 1. Therefore, where point D drops upon 
cutoff as the rarefaction wave tries to lower the absolute pressure at the gage location is 
depended on one of three levels [9]: 

E: If pressure at E is greater than vapor pressure 

F: If pressure at F is great than cavitation but less than vapor pressure 

G: If pressure at G is less than cavitation pressure 
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3. 


Bulk Cavitation 


The arrival of rarefacted or reflected wave near surface, also called cut-off due to 
water’s inability to withstand tensile force at air-water interface, creates BC near air- 
water interface. In general, unlike seawater, pure homogenous water with no impurities 
is able to withstand a considerable amount of tensile stress. Due to abundance of obvious 
biological and non-biological impurities, seawater can withstand very little tension. As a 
result, the layer of sea water at or near the water-air interface ruptures and the pressure 
falls to the vapor pressure of water, assumed to be approximately zero for all practical 
purpose when dealing with UNDEX [10]. Within the areas of rupturing water layer near 
the surface, a collective spallation takes places vertically and tensile stresses at its lower 
boundary tending the water particles to return to its original, pre-shocked position cease 
to exist. Of note, the motion of spallation is dictated primarily by the atmospheric, 
hydrostatic pressure and the acceleration due to the gravity. The radial and vertical 
thickness of typical BC zone for a 601bs charge of HBX-1 at a depth of 30ft is depicted in 
Figure 12. As shown, this charge size and depth creates a BC that is 574ft in radius and 
24.6ft deep at its maximum thickness; hence, as shown, its radius is several magnitudes 
longer than its radius. 



Figure 12. BC Zone for 60 lbs of HBX-1 at 30ft 


Next, the spalled water layer is driven back downward by force of gravity so that 
ultimately it must impact the underlying water and produces a traveling source of 
secondary pressure wave in the water called “hammer.” This traveling source and 
hammer created are almost equivalent to a “sonic boom” in the atmosphere and because 
of its focusing effects, the secondary pressure due to hammer can reach intensity as high 
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as initially propagating SW [10]. Interestingly, the region of cavitation between the 
upward traveling spalled water layer and the relatively quiescent underlying body of 
water have an average specific volume greater than that of seawater and this intervening 
region is termed the BC region due to its size and can be observed in figures 12 and 13. 



Figure 13. Photograph of Upper Surface of Bulk Cavitation, from [11] 


Finally, in an attempt to summarize numerous BC studies conducted and theories 
derived over the years since the 1960s, the critical sections within the technical paper 
published by Costanzo and Gordon in May 1983, titled “A Solution to the Axisymmetric 
Bulk Cavitation Problem” [9], as well as, sample case studies are all included in the 
Appendix of this study as a reference and also to be used in subsequent chapter’s analysis 
of LOD UNDEX phenomenon for its underlying characteristics of BCs. 

4. Motion of Gas Sphere and Its Empirical Bubble Formula 

During the expansion and collapse of cyclic bubble, the initial peak pressure in 
the expanding gas sphere is decreased dramatically due to energy dissipation to its 
surrounding from detonation. Flowever, the immediate state of the spherical bubble that 
follows still has much higher pressure than the equilibrium hydrostatic pressure at that 
depth [6]. This sphere or bubble has a tremendous outward velocity as its diameter 
increases rapidly and is proportional to the pressure magnitude at that exact time. The 
expansion of bubble continues for a relatively long period of time compared to other 
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transient UNDEX phenomenon like initial shock waves or BC. As this bubble reaches its 
maximum radius, the internal gas pressure decreases dramatically, while the expanding 
motion continues due to the inertia of water flowing outward [6]. As the internal gas 
pressure eventually falls below the equilibrium pressure (hydrostatic + atmospheric), this 
outward motion or flow stops and the bubble begins to contract at an increasing rate. The 
contraction or inward motion of the bubble continues until it reaches the maximum 
compressibility of the gas and reverses the motion once again due to reloading of 
pressure differences within the bubble and its surrounding. Hence, the created elastic 
properties of the gas and water due to detonating pressure and hydrostatic influences of 
surrounding water, in time creates a condition of oscillating system that repeats itself 
until it reaches the surface of the water [6], This cyclic nature of bubble expansion and 
collapse as function of time, as well as, its displacement as it travels up towards the 
surface due to hydrostatic forces induced by surrounding water and gravity is depicted in 
Figure 14. 



Figure 14. (a) Pulsation and (b) Displacement of Gas Bubble, from [6] 


The traveling cyclic bubble pulses and its comparison to corresponding pressure 
profile during UNDEX are shown in Figure 15. As mentioned, the initial shockwave 
pressure is the greatest at minimum bubble radius at the time of detonation. As the cyclic 
bubble pulse travels to near water-air interface, the detonating energies are dissipated and 
the bubble eventually reaches the surface to create plume into the atmosphere. This 
nature of cyclic bubble can induce hull whipping to surface vessels if period of the 
bubble pulse matches its lowest natural frequency; hence, the importance of analyzing 
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such phenomenon during UNDEX. Within the energy partition of the shockwave bubble, 
the initial shockwave dominates initial energy dissipation followed by cyclic bubble 
pulses (Figure 16). The initial SW energy can be as high as 50% of total bubble 
dissipated energy while the remaining energies are attributed to traveling cyclic bubbles 
pulses [8]. 



Figure 15. Traveling Cyclic Bubbles and Its Pressure Profile, from [5] 


Energy Partition of Shockwave Bubble 

■ Initial Shockwave 

■ 1st Pulse 

■ 2nd Pulse 


Figure 16. Energy Partition of Shockwave Bubble, after [8] 

Numerous theoretical methods and approaches to modeling and analysis of gas 
bubbles and its oscillations have been published over the years. However, for this study, 
a shallow water assumption was made to limit the scope of such analysis. As a result, 
two parameters of gas bubble emerge during UNDEX, its maximum radius and period of 
initial pulse [1]. Taking into account the correction factor that must be included due to 
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proximity of the air-water or water-bottom interfaces in shallow water, as well as, 
following the same intuitions in deriving shockwave pressure equations (2.5) and (2.6), 
the maximum bubble radius (A max ) in feet is derived [5]: 


A = K, 

max 6 


w 


D + 33 


(2.7) 


Next, the period of initial bubble pulse (T) is derived in seconds with its correction factor 
due to distance from surface, shown in Equation (2.8). The correction factor of 
numerical constant (a) is equal to 0.1 when Amax/D is less than 0.5 and when greater 
than 0.5, then a is approximated between 0.1 and 0.2 [8], 
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5. Secondary or Additional Pressure Pulses and Surface Effects 

As described previously, the bubble pulse due to cyclic expansion of gas as it 
travels to the surface of the water depends considerable on the water depth and proximity 
of the boundary surfaces at detonation. The first peak pressure of bubble pulse is usually 
no more than 10 to 20% of that of initial SW [6]. However, its duration is much longer 
and the areas under the pressure profile curve are considerably close [6]. During 
completion of each cyclic bubble pulse, a considerable amount of the energy initially 
present at detonation is lost due to energy transfer that occurs to its surrounding water as 
it travels and pulsate, as seen previously in figures 14 and 15. As a result, only the first 
bubble pulse is of practical significance since the successive pulses progressively weaken 
to negligent magnitude. The relations between initial SW and the first bubble pulse 
pressures and durations are shown in Figure 17 and Figure 18 shows close up details of 
the bubble pulse pressure generated from 300 pounds of TNT charges detonated at 
various depths in 100 feet total depth of water. 
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Figure 17. Pressure 60ft from 3001bs TNT Charge Detonated at 50ft, from [6] 



Figure 18. Various Bubble Pulses of 60ft at Various Detonated Depth, from [6] 

It is observed that the profile of the first bubble pulse curve in Figure 18 becomes 
increasingly irregular for initial charge position close to either surface or bottom [6]. 
Also, for shallow water or LOD conditions, the traveling cyclic bubbles and the reflecting 
pressure waves from the surface or bottom give rise to interference; therefore, as would 
be observed in this study, the later portions of observed pressure profiles are considerably 
different than what would be observed in an infinitely deep water domain [6]. Lastly, the 
interaction on the water-air surface boundaries as initial SW and bubble reaches the 
surface is also important. While the initial SW creates BC at this interface, the spray 
dome and subsequent bubble pulses can generate multiple overlapping plumes that 
literally shoot or hurl the mixture of explosive materials and water vapors into the 
atmosphere creating additional havoc during UNDEX. This surface phenomenon which 
marks the end of major UNDEX events is shown in Figure 19. 
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(a) Spray Dome (b) First Plume (cl) Second Plume 


Coused by 
Shock Wave 


Caused by 
First Bubble 
Pulse 


Caused by 
Second Bubble 
Pulse 


Figure 19. Surface Phenomena (a) Spray Dome or Spallation (Bulk Cavitation) 
Caused by Initial Shock Wave (b) First Plume Caused by the First 
Bubble Pulse (cl, c2) Subsequent Plume caused by Burping Bubble 

Pulses, from [5] 
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III. MODELING AND SIMULATION USING DYSMAS 


A. DYSMAS 

The software package used for this research is the DYSMAS finite element 
hydrocode. DYSMAS is a fully coupled; six degree of freedom (DOF) finite element 
(FE) based hydrocode that is designed to simulate three-dimensional (3D), UNDEX 
events to analyze fluid-structure interactions (FSA). The foundation of DYSMAS’ 
theoretical concepts is based on FE methods and like any common FE software suite 
currently available, it consists of two major fluid and structural solvers, Gemini and 
Dyna_N(3D), that is interfaced by the Standard Coupler Interface (SCI) along with a pre- 
and post- processor called, DYSMAS/P [12]. The basic DYSMAS suite architecture is 
depicted in Figure 20. 


DYSMAS Coupled Code Architecture 



Gemini 

Massively parallel explicit 
compressible flow solver 

Replay 

Rapidly replay structural loading 

Float 

Initialize draft, trim, & roll 

PreStress 

Statically load structure 



Dyna_N(3D) 

Massively parallel explicit 
finite element solver 


RigidBody 

Explicit rigid body dynamics 


Figure 20. Basic DYSMAS Coupled Code Architecture, from [12] 


The Gemini is an Eulerian fluid solver and the Dyna_N(3D) is a Lagrangian 
structural solver that evaluates the structural dynamic response during UNDEX. The SCI 
is the key coupler that connects and interfaces during a fully “coupled” simulation of both 
solvers by passing information between the two solvers at the end of each Eulerian time 
step, maintaining fully coupled interface. There are several other structural solvers that 
can be used in coupled DYSMAS run, such as RigidBody solver. DYSMAS/P and other 
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post processors like MATLAB were used for detailed analysis on the results, both 
graphically and numerically. 

1. GEMINI 

For this study, only the Gemini portion of DYSMAS was utilized in Eulerian fluid 
behavior analysis during UNDEX. The Gemini code is a finite-difference, Eulerian fluid 
equation solver that solves fluid mesh by running a higher order Godunov method 
algorithm to solve the fluid domain at each time step [13]. Designed specifically to 
simulate the shockwaves, cavitation and bubble jetting phenomena of UNDEX, Gemini is 
capable of solving flow fields involving several Eulerian fluid types. Using conservation 
of mass, momentum, and energy, the Godunov method algorithm is applied to solve each 
Eulerian equation for each cell within the fluid domain with a major assumption of Euler 
equation being that the flow is frictionless or inviscid. 



Figure 21. Basic Components of GEMINI, the Eulerian Solver, from [12] 


Gemini consists of several additional components codes: GemGrid, PreGemini, 
Float, and Prestress are preprocessors. The user defined and required “input” files 
(shown in blue, Figure 21) are used to specify each component codes. The sample input 
files are also included in the Appendix of this report. These sub component codes are 
used to customize fluids domain set up for a given UNDEX problem. GemGrid is used 
to set up the initial Euler cell grid, up to three-dimensions and can be used for any grid 
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refinement that is needed around specific areas of interest in the flow field, specifically 
the fluid bubble, cavitation and the Lagrangian structures. PreGemini is used to fill the 
Euler cell grid with the user-defined Eulerian materials that utilizes the given “material 
files” that encompasses specific material’s equation of state and initial properties, such as 
energy, density, and etc. In order to accurately position and prestress the Dyna_N 
structure in the fluid domain prior to transient analysis, the Float and Prestress processors 
are used. 

As Gemini solves each Eulerian cell grid, it dumps and saves all flow field and 
historical data into a bin file, which in turn are finally processed by GemHis and 
Gemfield processor. GemHis processes the user specified recorded UNDEX variables 
and locations within the fluid domain for analysis. Gemfield generates contour plots for 
the similar user defined UNDEX variables of designated plane snap shot at specified time 
increments. Once flow field and historical data files are generated using GemHis and 
Gemfield, DYSMAS/P, an independent DYSMAS pre/post processing suite or MATLAB 
can be utilized to further analyze the Eulerian and Lagrangian solutions. Samples of all 
input files are included in the Appendix of this report. 

B. DYSMAS/P AND MATLAB 

DYSMAS/P and MATLAB were used for pre and post processing of DYSMAS 
simulations. The preprocessing of all Lagrangian FE models can be conducted using the 
DYSMAS Pre-Processor 2010 and some Eulerian characteristics like BCs were done 
using MATLAB. The preprocessing capabilities of DYSMAS/P can be used for the 
creation of the FE model’s structure, assignment of material properties, boundary 
conditions, body forces, motions, and fluid-structure interface segment definitions. The 
structural model is then written into the specific input cards required Dyna N, in turn 
used as additional input files during coupled runs. The DYSMAS/P postprocessor allows 
the user to analyze the entire fluid-structure dynamic response, allowing user to extract 
specified data and a snap shot of specified time steps for graphical and numerical 
analysis. MATLAB was also used in numerical post processing of DYSMAS solutions. 
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C. BASIC DYSMAS SIMULATIONS SETUP 

The basic DYSMAS simulation set up entailed nearly the same basic steps and 
approaches on all simulated cases with minor differences depending on the type of run as 
fluid only or fluid-structure coupled run. For this study, fluid only runs were conducted. 
Following steps, in the order as listed and corresponding input files were utilized: 

1. Utilize DYSMAS/P and MATLAB for preprocessing as needed for 
Lagrangian FE structure model and specific Eulerian fluids characteristics 

2. Generate Euler grid using Gemgrid 

3. Fill the fluid domain using Pregemini & Material Files 

4. Check fluid field prior to Gemini execution using Gemfield 

5. Run Gemini (For coupled run, Gemini will be specified as such, invoking 
fully coupled 3D run using Dyna_N(3D)) 

6. Process Gemini/Dyna_N(3D) solutions using Gemhis and Gemfield 

7. Detailed post processing and qualitatively validate results prior to using 
DYSMAS/P and MATLAB for final post processing 

Like any other FE methods or software, the “art” of utilizing FE was also 
considered throughout the process, especially during the mesh or grid set up, validation 
of initial results and distinguishing the best use of computational power versus time for 
each simulations. 

D. EULER GRID, FLUID DOMAIN SETUP AND CHARGE PARAMETERS 

During grid setup using Gemgrid, the most important decision made were where 
to refine grid and where to maintain as default coarse grid of approximately 0.656 ft per 
grid cell. According to Didoszak and from previous studies conducted at NPS, 
DYSMAS simulations are highly sensitive to the mesh size surrounding the explosion 
[14]. It was found that the smaller the mesh or grid size, the closer in accuracy in the 
simulation for both empirical and experimental predictions for various UNDEX 
parameters. Obviously, these assumptions had to be taken into consideration with a grain 
of salt, because mesh refinement is not the solution to all FE problems. In addition, 
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Gemini’s free boundary conditions did not adequately allow fluid to flow in and out of 
fluid domain quickly enough, if the fluid domain was too small in comparison to the 
maximum bubble diameter. Lastly, when setting bottom boundary conditions, anything 
less than “wall” or rigid gave unexpected errors in flow field and historical results. This 
finding and realization was particularly useful and had to be considered in shallow water 
simulations and solutions where overall fluid depth was restricted considerably compared 
to deep open ocean scenarios. 


/ Air Column 

-» x 

> Charge (X, Y, Z) 

Water Column 


_ Clay Column 

Figure 22. Typical Fluid Domain Setup for Gemini 

Following the grid set up, the fluid domain was initialized and populated through 

the Pregemini input file. As typical Eulerian fluid domain setup is depicted in Figure 22 

and the critical step in this procedure was the seeding of the fluid domain with the correct 

materials at the appropriate states and depth or geometry and verifying it with Gemfield 

prior to Gemini execution. In order to populate the fluid domain with different materials 

such as explosives, air, water, etc. Pregemini was utilized with given various material 

files. Within the material files contained specified equation of state and variables 

(pressure, energy, and density) that can also be manipulated, if needed, for specified 

materials at user’s discretion. As part of the fluid domain set up within the Pregemini 

input file, except where specifically noted, a standard charge of HBX-1 was used. The 
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charge weight and depth were varied or kept constant according to the specificity of each 
simulation conducted. 

E. SUMMARY OF SIMULATION SETUP 

Taking all previously discussed factors of UNDEX theory and lesson learned 
from initial simulations and previous studies, the deciding values, properties and setup for 
each simulation started from the basic guidance and logics established in this section. 
This basic set up was then further explored and developed into complex and more 
realistic Eulerian fluid model as each DYSMAS solutions were qualitative validated in 
accordance with UNDEX theories discussed in previous sections of this report prior to 
final simulation runs. Each case studies and simulation set up, including its parameters 
and contributing factors can be found throughout subsequent chapters of this report and 
the complete outline of case studies can be found in the Appendix of this report. 
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IV. INTRODUCTION LITTORAL OCEAN DOMAIN 


Recent proliferation of asymmetric threats, such as, diesel subs, mines, fast patrol 
boats, coastal integrated defense in littoral areas and even drones makes littoral ocean 
domain (LOD) a battle space that USN must carefully consider and analyze in the future 
maritime battle scenarios [15]. The term, “littoral” water (as depicted in Figure 23) has 
been used to extensively to describe such domain, where shallow water and its depth, as 
well as, shoreline and possible fresh water composition and tidal waves play a critical 
role in battle space management and the type of vessels or forces that USN can deploy to 
this domain. 



Figure 23. Littoral Ocean Domain, from [16] 


The term “Littoral” and “Shallow” water domains are synonymous for this study 
and any body of water not of open-ocean, less than or equal to 500 feet depth and close to 
shore is in this category. Figure 24 shows the vast majority of ocean domain as open 
ocean and from a depth perspective, the LOD is only a fraction of water domain above 
the continental shelf. However, the vast majority of economic transaction and military 
conflicts are conducted in LOD and this number continues to grow [17]. Taking it 
further, the LOD will be bounded generally anywhere but open ocean and within the 
domains of depth ranges from 150 to 300 feet or approximately 10 to 20 times the 
average draft of USN’s Littoral Combat Ship (LCS, 12.8 feet / 3.9 meters) and presence 
of tidal currents, gradual decent of ocean bottom shelf to open ocean depth, presence of 


29 










sedimentary debris throughout water and numerous bottom sediment types to include 
clay, sand and rocks are some of many characteristics that LOD embodies. 



Figure 24. World Ocean Domain, from [18] 


A. LITTORAL OCEAN DOMAIN RELATED TO UNDEX 

To further simplify LOD UNDEX problems, only two types of bottom sediments 
will be considered for this study—clay and sand—for its reflectivity or absorption of 
energy when induced to forces incurred during UNDEX. The distinct differences in 
behavior of cohesive clay bottom type and non-cohesive sand will not be explored during 
this study, rather its reflectivity or absorption of energy. Sand in general will absorb 
more energy from UNDEX due to its increased porosity between sand particles. Also, 
the angle of decline or incline in shallow water areas, like beaches or shoreline coming in 
from the deep abyss will be assumed to be small or negligible. Therefore, only the depth 
and side boundary conditions will be explored for UNDEX phenomenon in this research. 

There are four primary forces that can be induced to a surface vessel operating in 
an LOD. They are initial detonating shock pressure wave, the “hammer” shockwave 
created by the BC formation and dissipation, the bottom bounce (BB) or laterally 
reflected (LR) pressure waves and rising cyclic gas bubbles due to explosive materials. 
While all four forces are similar to that of what is experienced in an open ocean, the 
effect of bottom bounce, LR pressure wave and additional shockwave due to cavitation in 
LOD is what sets this domain apart. The rise and effect of gas bubble is purely 
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determined by the size of initial bubble radius; therefore, depending on the charge size 
and initial bubble radius, the subsequent cyclic bubbles may or may not exist. The most 
distinctive feature that differs in LOD BC formation or dissipation is the follow on 
disorderly, additional expansion and collapse of BC. For extremely shallow depth below 
100ft, the bottom-reflected shock or BB arrives before cavitation collapse starts, and as a 
result, it interacts with BC in a complex manner, narrowing its region and hastening its 
collapse to its original state [11]. Additionally, it induces additional cavitation in the 
vicinity of the bubble and gives rise to a more disorderly collapse involving the original 
cavitation region, as well as, additionally induced cavitated areas. The shallower depths 
also play a role in formation of and subsequent cyclic behavior in the detonating bubble. 
Again, if the depth of LOD is smaller than the radius of first period cyclic wave of 
detonation, the debris within the detonating sphere will spew out in its first expansion 
creating a one large plume rather than several subsequent plumes. 

B. UNDEX BENCHMARK PROBLEM FOR LITTORAL OCEAN DOMAIN 

Using DYSMAS, a benchmark UNDEX problem for this study is created and 
validated (Figure 25). This fluid or Eulerian domain contains three distinct layers of 
Eulerian materials to include clay bottom, sea water and a layer of air. The origin of axis 
lies at water-air interface right above the charge, and the depth and lateral distances of 
domain are varied in order to meet specific objectives for each run cases. Both charge 
and target locations are prescribed at some point in the domain, (X, Y, Z). In addition, 
HBX-1 will be used as the charge for the remainder of this study, due to carefully 
documented characteristics of its composition and UNDEX behaviors. HBX-1, also 
known as High Blast Explosive was developed during WWII as desensitized 
modification of Torpex explosives and its composition includes 38% TNT, 40% RDX, 
Aluminum Wax and Lecithin [19]. RDX, also known as Research Department Explosive 
is an explosive nitroamine widely used in military and industrial applications. 

The setting of boundary conditions for each run was of critical importance; since, 
they detennined the approximated interactions between various interfaces of UNDEX 
geometry and the way shockwaves and energy will travel accordingly. The ocean bottom 
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and lateral boundary (B1 though B4), such as channel walls, were assumed to be a rigid 
body initially to simplify the modeling requirements. For DYSMAS, there are essentially 
three boundary settings: free, wall and 0 to 99% reflectivity. The condition “free” 
includes gravity correction, relevant for the z-direction only. The condition “wall” 
specifies a reflecting or rigid boundary condition. Lastly, the numbers between 0 and 1 
can be used for partially reflecting boundary conditions [12]. 



* 


B3 


Figure 25. UNDEX Benchmark Problem, Defining Littoral Ocean Domain 

To confirm initial assumptions made regarding LOD in reference to UNDEX 
parameters, several initial case studies were conducted with DYSMAS by first simply 
varying ocean depths. Initial depths were varied from 1000ft to 100ft, and additional 
depth ranges were further broken down and simulated to observe or validate UNDEX 
characteristics within that range. The simulated cases from 1 through 14 are summarized 
below in Table 1. The UNDEX characteristics were analyzed using several results 
obtained from each run to include: flow field, shockwave characteristics, vertical take¬ 
off velocity and BC. 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,Z) [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[«] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

1 

2D, Euler 

800,1,1000 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

2 

2D, Euler 

800, 1, 900 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

3 

2D, Euler 

800, 1, 800 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

4 

2D, Euler 

800, 1, 700 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

5 

2D, Euler 

800,1, 600 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

6 

2D, Euler 

800,1, 500 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

7 

2D, Euler 

800, 1, 400 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

8 

2D, Euler 

800,1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

9 

2D, Euler 

800, 1, 200 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

10 

2D, Euler 

800,1,175 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

11 

2D, Euler 

800, 1,150 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

12 

2D, Euler 

800,1,125 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

13 

2D, Euler 

800, 1, 100 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

14 

2D, Euler 

800,1, 75 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 


Table 1. Case Studies for Defining Littoral Ocean Domain 


1. Flow Field Analysis 

The overall flow field results for deep (greater than 500ft) and shallow (less than 
300ft) ocean, or LOD, are depicted in figures 26 and 27. Both flow fields depict pressure 
contour plots or “flow field” at time steps of (a) 10msec (b) 35msec (c) 50msec and (d) 
75msec. These time steps were chosen to coincide with the most significant UNDEX 
events taking place. Both figures 26 and 27 at time (a) 10msec show identical initial 
shockwave expansion as it intercepts the water-air interface. Once in contact, the 
characteristic of both cases take a different turn. The expansion and initiation BC can be 
clearly seen in subsequent field plots at time 35 and 50msec, with initial nucleation of 
collapse (Figure 26(b)) and hammer shockwave (Figure 26(c)). The traveling hammer 
shockwave can also be seen in Figure 26(d) as it propagates deep into the ocean. At this 
point, the same hammer shockwave have already reached the surface and dissipated. As 
mentioned in previous chapter, hammer shockwave is magnitude of the pressure pulse 
generated by the collision of the accreting lower and upper cavitation boundary at 
closure. 
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Figure 26. 500ft Ocean Depth at Time Step (a) 10 (b) 35 (c) 50 (d) 75msec 


The flow field results for 200ft ocean depth are depicted in Figure 27. There are 
several contrasting differences to previous 500ft depth flow field results. First, the 
bottom bounced (BB) reflected shockwave are shown in Figure 27(b). The magnitude of 
this BB shockwave at reflection is identical to initially radiating shockwave. As it travels 
back up towards the surface, the BB shockwave intercepts the “hammer” shockwave 
from collapsing BC, shown in Figure 27(c). Last, as BB shockwave passes through the 
hammer shockwave and hit the water-air interface, it invokes additional BC that can be 
seen in Figure 27(d). Clearly the chaotic nature of shallower depth UNDEX can be 
contrasted from the 500ft flow field results and start to emerge. 
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Again, the 500ft depth UNDEX results were representative of all depth cases 
deeper than 300ft and 200ft depth results were representative of cases shallower than 
300ft. From these flow field results, there are four differing characteristics of deep and 
shallow water UNDEX. First, the way in which initial shockwave travels and interacts 
with “other” reflecting waves differ. Unlike open ocean, as initial shockwave reflects off 
shallow bottom, LB or water-air interface and heads towards a point of intersection of all 
three waves, additional cavitation and shockwaves are created. Second, the BC collapse 
of shallow water is more violent or higher in pressure fluctuations during “hammer” 
shockwave propagation due to additionally induced cavitation or shockwaves. Third, the 
BB shockwave during shallow water UNDEX induces additional energy back into the 
fluid field creating additional BCs and shockwaves that can be as close as the initial 
shockwave and energy. The additional BC can be seen in Figure 27(d) as the BB 
shockwave crashes into the water-air interface, inducing another set of spallation and 
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additional formation of cavitations. Lastly, the additional BCs that are formed during 
BB’s collision of water-air interface create an additional “hammer” shockwave that also 
can be as high as the initial closure pulse. Another flow field results are shown in Figure 
28 for shallow water UNDEX to show the violent and multiple cavitations forming due to 
interactions of several reflected shockwaves from the bottom. As shown, the white 
portion of flow field indicates BC formations. Clearly, the additionally induced BC 
formation and its chaotic behavior can be seen and at times greater than that of initial BC 
in size and scope. 



c. t s ♦ 51.58 ms d. t s + 121.6 ms 

Figure 28. UNDEX of 11021bs TNT at Charge Depth of 66ft, Water Depth of 

131ft, from [11] 


2. Target Pressure Analysis 

The open ocean and LOD’s target pressure profile during UNDEX are depicted in 

figures 29 and 30. Both figures show that the initial shockwave pressure at target 

location to be approximately 55psi. In Figure 29, this initial shockwave pressure is 

followed by an immediate cutoff pressure at 8psi with reverberation within the layer of 

target depth and air-water interface, and finally a second shockwave peak pressure at 

approximately 22psi representing the “hammer” shockwave pressure during the initial 

BC closure. To get a perspective on shockwave pressures, 14.7 psi is the atmospheric 

pressure that is felt at sea level and 55 psi of initial peak pressure is from a 251bs charge 

of HBX-1 at a depth of 50ft of sea water. 55psi is close to what you will feel at 
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approximately 120ft deep water, stationary. The rest of the time period show negligible 
pressure peaks compared to initial and hammer shockwave, although are clearly present. 
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Figure 29. Target Pressure Profile (300 to 1000ft Ocean Depth) 


As seen on the overall flow field plots previously, Figure 30 shows additional 
pressure waves generated due to BB waves and additional BC closures. For 75, 100 and 
125ft of water depth, the additional shockwave present due to BB and additional BC are 
found to be up to 90% of what initial shockwave experienced. For the 100ft case, the BC 
collapse pulse and BB SW that follow are almost identical. The pressure fluctuations are 
clearly more observable throughout the course of the shallow pressure time domain. The 
pressure fluctuations felt as individuals may or may not cause structural damage 
depending on the magnitude and material properties of structures present. However, 
accumulation of dynamic nature of additional shockwave can and will increase the 
likelihood of the structural damage due to dynamic loads. 
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Target Pressure Profile (75 ~ 200ft Ocean Depth) 



Figure 30. Target Pressure Profile (75 to 200ft Ocean Depth) 

3. Vertical Take-Off Velocity Analysis 

Next, the vertical take-off velocities (VTO) are analyzed for the open ocean and 
shallow water UNDEX events (figures 31 and 32). In short, VTO represents magnitude 
of z-direction velocity experienced at target location during UNDEX. The abrupt 
increase in velocity at a constant rate followed by a varying decay of magnitude shows 
close resemblance to initial shockwave characteristics. The VTO for ocean depths of 300 
to 1000ft is shown in Figure 31 while Figure 32 depicts VTO for ocean depth of 
shallower depths, less than 300ft. Both VTO results show the peak of velocity at 1.75ft/s 
before decaying back close to its stationary state or into negative velocity due to 
additional fluid movements created by additionally reflected shockwaves or formed 
cavitations. For deeper water, rest of the remaining VTO during period of UNDEX 
seems to decrease, approaching close to negative 2ft/s. 
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Vertical Take-Off (VTO) Velocity (300 ~ 1000ft Ocean Depth) 



Figure 31. Vertical Take-Off Velocity (75 to 200ft Ocean Depth) 


An entirely different phenomenon takes place for shallower depth in Figure 32. 
Although the peak VTO is identical as mentioned, the additional VTOs felt due to BB 
and additional BC formation and collapse can be seen. Again, closely resembling 
shockwave pressure, the addition VTO experienced at target locations is found to be up 
to and even as greater than 50% of the initial VTO. The numerous additional VTO peaks 
and chaotic nature of shallower depth is apparent and closely resembles previously 
analyzed flow fields and shockwaves. In short, shallower the depth, VTO is felt more 
violently and duration is longer. As observed in pressure profiles previously, the distinct 
characteristics of less than 300ft starts to emerge for both UNDEX parameters. 
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Vertical Take-Off (VTO) Velocity (75 ~ 200ft Ocean Depth) 



Figure 32. Vertical Take-Off Velocity (75 to 200ft Ocean Depth) 


4 . Bulk Cavitation Analysis 

Lastly, the BC’s overall characteristics and size are analyzed. The total volume of 
BC formed during open ocean and shallow water UNDEX are shown in figures 33 
through 35. The open ocean cases, Figure 33, show a peak BC volume of approximately 
6x10 ft (4.5 x 10 Gallons). This is a tremendous amount of water moved by a 251bs 
charge of FIBX-1; hence, like the initial shockwave of 55 psi, one can gauge at the 
tremendous energy that is imparted to the water during just the initial part of an UNDEX 
event. For all 300 to 1000ft depth oceans, or cases without bottom reflection, the peak 
and total BC volumes are almost identical. By 35msec, BC starts to collapse and by 50 to 
60msec, complete dissipation of initial BC can be observed both in flow field plots 
previously and BC volume plots. 
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Bulk Cavitation Volume (300 ~ 1000ft Ocean Depth) 



Figure 33. Bulk Cavitation Volume (300 tolOOOft Ocean Depth) 


The shallow water BC characteristics paint a different story as briefly mentioned 
in previous section. The increased cavitated zones induced by the BB wave and are 
clearly evident in shallow water (figures 34 and 35). From 75 to 200ft of water, the total 
BC formed in 75ft of water surpasses all other depth’s BC formed individually, as shown 
with red dotted horizontal line in Figure 34. In 100ft of water, the BC peak is observed to 
be 6 x 10 ft , over four million gallons. Four million gallons is approximately equivalent 
size of present day water tower that exists throughout various cities. This dramatic 
increase in additional BC formation for shallower depth can be attributed to many factors, 
but for this study, focus will be on characterizing shallow depth and effects of boundaries 
at the bottom and lateral position of water domains that are present in shallow ocean or 
LOD. 
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x IQ 5 Bulk Cavitation Volume (75 ~ 200ft Ocean Depth) 



Figure 34. Bulk Cavitation Volume (75 to 200ft Ocean Depth) 


Bulk Cavitation Volume (100 ~ 200ft Ocean Depth) 



Figure 35. Bulk Cavitation Volume (100 to 200ft Ocean Depth Cases) 
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V. INVESTIGATION OF CHARGE SIZE AND DEPTH 


Before moving further, the effects of charge size and depth must be addressed and 
validated for LOD. To do so, the same UNDEX fluid domain was used as in previous 
cases, except the ocean depth was maintained at 500ft deep while charge sizes and depths 
were varied. The depth at 500ft was chosen to analyze ocean depth of greater than 300ft, 
but not quite deep enough to see the negligible changes in UNDEX characteristics caused 
by charge size and depth. The ocean domain for this case studies are shown in Figure 
36 and Table 2, where cases 15 through 24 are summarized for its detailed set up. 


Z 
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B3 


Figure 36. Littoral Ocean Domain for Charge Size and Depth Case Studies 


Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,Z) [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

15 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

16 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

200 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

17 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

300 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

18 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

400 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

19 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

500 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

20 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

100 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

21 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

200 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

22 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

300 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

23 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

400 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

24 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

500 

Clay 

Free 

Free 

Wall 

Wall 

0.656 


Table 2. Case Studies for Varying Charge Size and Depth 
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A. 


VARYING CHARGE SIZE (100 - 500LBS OF HBX-1) 


1. Flow Field Analysis 

The overall flow field result of varying charge sizes are shown in Figure 37. The 
time step chosen represents three distinctly important events taking places: (a) 
propagation of initial shockwave (b) maximum radius of BC and (c) BC collapse where 
hammer shockwave is generated. Again, as the initial shockwave slams into the water-air 
interface, the BC develops, maturing into maximum radius by 55msec and eventually 
collapse at 75msec as the initial shockwave continues to propagate into the depth and 
remaining ocean domain. The chosen 500ft ocean depth does not interfere with initiation 
and collapse of BC and secondary and tertiary UNDEX events are not readily observed 
during the remainder of chosen time step of 150msec, hence not included. For rest of the 
charge size cases, the differences in the flow field results were maximum BC radius and 
time of BC closure, obviously due to its sizes. 
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Figure 37. Flow Field at Time Step (a) 10 (b) 55 (c) 75msec 


2. Target Pressure Analysis 

The results of target pressure for varying charge sizes are depicted in figures 38 
and 39. Figure 38 shows the maximum pressure of 5001bs FIBX-1 at approximately 
l,000psi. Next, Figure 39 shows maximum pressure peaks of 100 to 4001bs HBX-1. The 
difference of peak pressures between 400 and 5001bs HBX-1 is close to 700psi, while 
differences of peak pressure between 100 and 4001bs HBX-1 is close to 200psi. These 
drastic and almost exponential differences in pressure peak between 500 and 4001bs 
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charge is of significant, since it can be used to bound problems via similitude to pick an 
ideal weight that is most suitable for a given case studies. 


Target Pressure for Various Charge Sizes 
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Figure 38. Target Pressure Profile (100 to 5001bs Charge Weight) 


Target Pressure for Various Charge Sizes 



Figure 39. Initial Target Pressure Profile (100 to 5001bs Charge Weight) 
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3. Vertical Take-Off Velocity Analysis 

The VTO results for various charge sizes are depicted in Figure 40. The 
maximum VTO for 5001bs peaks at approximately 3 7ft/s while 4001bs or less charges 
peak at 12ft/s or less. Similar to pressure peaks observed previously, the differences in 
the VTO between 100 to 4001bs charges and 5001bs charges are quite drastic. Other 
characteristics that are noticeable are minor VTO that starts to take shape after 50msec 
for 4001bs or less charge sizes. These minor VTO can be attributed to the smaller BC 
collapsing and inducing additional fluid movements that can be observed during the time 
step of simulations. For 5001bs charge, the sheer size fluids moved and BC formed 
during detonation takes a lot longer to return to its original state than the designated 
150msec time frame, hence, the magnitude of its VTO. 


Vertical Take-Off Velocity (VTO) for Various Charge Sizes 



Figure 40. Vertical Take-Off Velocity (TOO to 5001bs Charge Weight) 


4 . Bulk Cavitation Analysis 

The BC formation and dissipation for various charge weight is shown in Figure 
41. Again, the drastic differences in BC volume between 100 to 4001bs charge weight 
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compared to 5001bs charge cannot be ignored. The shape and speed of cavitations are 
consistent with previous cases studied. The sheer amount of water that is cavitated due to 
energy imparted into the water via propagating initial shockwave is close to 3x10 ft 
(2xl0 8 Gallons) for 5001bs charge, while 4001bs charge cavitates close to 7.5xl0 6 ft 3 
(6xl0 7 Gallons) water domain. It is important to keep in mind that at this point, the 
effect of BB or LB’s reflected waves and their cause of additional cavitation are not even 
considered due to chosen 500ft ocean depth. At this depth, for the chosen charge weight, 
the BC collapses prior to any interfering additional reflected shockwave present. 
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Bulk Caviation Volume for Various Charge Sizes 
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Figure 41. Bulk Cavitation Volume (100 to 5001bs Charge Weight) 


B. VARYING CHARGE DEPTH (100 - 500FT OF 100LBS HBX-1) 

1. Flow Field Analysis 

The overall flow field result of varying charge depths are shown in figures 42 and 
43. Figure 42 depicts a 300ft charge depth case and Figure 43 depicts 500ft charge depth 
case where the charge itself is at the water-clay (bottom) interface. The time step chosen 
to represent four distinctly important events taking places are: (a) propagation of initial 
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shockwave (b) initial shockwave reaching water-air interface (c) maximum BC radius 
and (d) additional cavitation formed due to BB shockwaves. 



Figure 42. Flow field at Time Step (a) 10 (b) 60 (c) 73 (d) 150msec (300ft 

Charge Depth) 

Unlike previous cases, the shockwave in Figure 42 takes a lot longer to reach the 
water-air interface due to its charge depth. As a result, the overall pressure, VTO and BC 
events are delayed up to 60msec. Further, the shockwave pressure and its magnitude that 
reached the water-air interface has drastically reduced, shown in Figure 42(b), since its 
energy has been continuously dissipating into the water as it travels up to the surface. 
The usual BB shockwaves that are reflected passes the charge depth (Figure 42(c)) and 
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reaches water-air interface after the BC have collapsed. The collapse of additional 
cavitation formed by BB shockwave can be seen in Figure 42(d). Next, for 500ft charge 
depth, two distinctly different events take place, shown in Figure 43. First, since the 
charge is embedded into the bottom clay, as detonation takes place, the initial shockwave 
and BB shockwaves are merged to form a single shockwave that propagates up to the 
surface of the water. As this emerging shockwave travels up to the water-air interface 
and collides, the BC fonned is much greater than any other charge depth combined. 



Figure 43. Flow field at Time Step (a) 10 (b) 102 (c) 119 and (d) 130msec 

(500ft Charge Depth) 
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From the flow field results, this huge BC can be seen in Figure 43(c) and (d). As 
this BC continues to grow, a similar characteristic can be observed in terms of additional 
cavitation fonned during shallow water cases observed previously. In short, due to a 
combined shockwave in lieu of two initial and BB shockwaves, the BC formed for this 
case is almost double in its original size if the charge depth was deeper than its half-way 
point of water domain or submerged into the bottom clay. For charge depth closer to 
water-air interface, similar phenomenon can be observed, except the absence of BC due 
to 50% of detonation energy expanded directly into the atmosphere. 

2. Target Pressure Analysis 


Target Pressure for Various Charge Depth 



Figure 44. Target Pressure Profile (100 to 500ft Charge Depth) 


The overall target pressure profiles are depicted in Figure 44 for various charge 
depth cases. The pressure peak of 100ft charge depth shows maximum pressure of 
approximately 135psi and 500ft charge depth shows a maximum pressure peak at 
approximately 40psi. As mentioned previously in overall flow field results, as the charge 
depth is increased and with it, the distance to water-air interface is also increased; hence, 
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the actual magnitude of shockwave’s energy that collides with the interface is greatly 
reduced since it must travel through greater volume of water. As a result, the total 
volume and timing of BC formation is greatly reduced and delayed, except for 500ft 
charge depth case. Taking closer look at the 500ft charge depth case, the secondary 
shockwave is much greater than that of its counterpart (Figure 45). As described 
previously through overall flow field results, this is due to the emergence of two 
shockwaves: initial shockwave and BB shockwaves. This secondary shockwave is 
similar in magnitude of 100ft charge depth cases. The tertiary and rest of the shockwaves 
observed are consistent in magnitude and time duration except for 500ft charge depth 
case. For this case, tertiary and rest of the shockwaves are actually missing since these 
shockwaves were combined to form greater secondary shockwave. 


Target Pressure for Various Charge Depth 



Time [msec] 


Figure 45. Close-Up Target Pressure Profile (100 to 500ft Charge Depth) 


3. Vertical Take-Off Velocity Analysis 

The VTO for various charge depth cases are depicted in Figure 46. Consistent 
with previous UNDEX parameters, the decrease in VTO as charge depth is increased can 
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be observed. The 100ft charge depth’s VTO is close to 5ft/s while 500ft charge depth is 
shown to exhibit VTO of close to 1.75ft/s. Not surprisingly, the differences to the 400ft 
charge depth’s VTO is minimal when compared back to the bottom depth case of 500ft. 
Another characteristic of VTO that can be observed is that as the initial VTO comes back 
to the original state, this overall state is increased for shallower depth. For example, 
VTO for 100ft returns back to greater than 0.5ft/s and hovers at this state of velocity 
through the remainder of its timeframe. This characteristic is also true for the remaining 
cases, except for the bottom depth case, as pointed out in Figure 44. Instead of initial 
VTO spike and returning to constant secondary state, there is a secondary spike in VTO 
followed by tertiary fluctuated spikes. Again, these secondary and tertiary spikes are 
attributed to the additional cavitation formation at initial shockwave interception of 
water-air interface, rather than later time for rest of the depth cases. Lastly, the tertiary 
VTO spike settles out at a lower magnitude than what the previous cases exhibited. 


Vertical Take-Off Velocity (VTO) for Various Charge Depth 



20 40 60 80 100 120 140 160 180 200 

Time [msec] 


Figure 46. Vertical Take-Off Velocity (100 to 500ft Charge Depth) 
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4. Bulk Cavitation Analysis 

The BC volume formation is depicted in Figure 47. The overall fonnation and 
dissipation of BC is consistent with what was observed in previous cases where BC is 
initiated at initial shockwave’s interception of water-air interface, reaches its maximum 
radius at some point depending on the charge size and additional cavitations are formed 
depending on the presence or lack of bottom. The BC volume for charge depth of 100 
and 200ft follows this consistent pattern observed previously; however, the rest of the 
cases differ. For a charge depth of 300 and 400ft, an additional “hump” can be observed 
at the posterior section of BC formation and dissipation. This is due to BB shockwave 
catching up to the BC dissipation and as it interferes with the BC collapse, the BB 
shockwave induces additional cavitation to be formed. This phenomenon is greatly 
amplified in the 500ft charge depth case where initial BC and additional cavitation are 
formed simultaneously, causing greater volume of overall BC. As observed, the overall 
max BC formed for the bottom depth case is close to 2.3x10 ft (1.7x10 Gallons), while 
the 100ft charge depth case’s overall max BC volume is drastically less than the fonner at 
1.8xl0 6 ft 3 (1.3xl0 7 Gallons). While the amount of overall BC formed is not as 
important as additional cavitations formed due to additional shockwaves that can be 
created, the sheer volume of additional BC created due to charge depth differences cannot 
be ignored. Hence, the collapse of BC formed by the 500ft charge depth will have 
greater residual shockwaves originated from its hammer pulse. 
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Bulk Caviation Volume for Various Charge Depth 



Figure 47. Bulk Cavitation Volume (100 to 500ft Charge Depth) 
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VI. INVESTIGATION OF LATERAL BOUNDARY CONDITIONS 


In LOD, the lateral boundary (LB) condition (B2) play a role in truncating BC 
formation, as well as, drastically enhancing the role of reflected wave, much like the BB 
effect observed previously. Cases such as channels, test ponds and areas of LOD where 
lateral spaces are limited, are critical to in-depth analysis of this case due to its effects on 
UNDEX parameters. To take a look at such phenomenon, the placement and reflectivity 
of lateral wall boundaries were varied during UNDEX (Figure 48). The placement of LB 
conditions (B2.1, 2.2, and 2.3) were a matter of where the BC’s edge was located and 
accordingly placed inside, at the edge or outside of the maximum BC radius. These 
represent the physical boundaries that the target and charge placements are as prescribed 
in the Case 25 to 28 in Table 3. In addition, reflectivity of boundary case 26 (B2.1), 
within the maximum BC radius, was varied to analyze its effect and influences on 
UNDEX parameters. 
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Figure 48. Lateral Wall Boundary Conditions (B2.1, B2.2 and B2.3) 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 

(X,Y,Z) [ft] 

Target 

Depth 

[«] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

25 

2D, Euler 

230, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

26 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

26.1 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

0.50 

Wall 

Wall 

0.656 

26.2 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

27 

2D, Euler 

197, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

28 

2D, Euler 

230, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 


Table 3. Case Studies for Varying Lateral Boundary 


A. VARYING LOCATIONS 

The next set of numerical results show initial and follow on shockwave pressure, 
vertical take-off velocity and BC volume experienced at target location and depth. As 
mentioned in Table 3 for case studies 25 through 28, the lateral target location is at 100ft 
and 0.820ft deep, for all practical purpose in regards to mesh size at 0.656ft, this depth is 
at shallow surface depth near water-air interface. 

1. Flow Field Analysis 

The overall pressure flow field results are shown in figures 49 to 52. Figure 49 is 
the base case where LB was set well outside the maximum BC radius with free boundary 
condition. The snapshots of flow field were taken at time step 10, 30, 50 and 111msec 
where significant UNDEX events took place. Consistent with previous cases, the flow 
field results show expected behaviors of (a) radiating initial shockwave, (b) initiation and 
dissipation of BC, (c) BC’s “hammer” shockwaves and (d) influence of other reflected 
waves from bottom or lateral boundaries. In Figure 49(d), the BB SWs are shown 
passing the charge depth as it intercepts the water-air interfaces. The “free” LB condition 
set at x-maximum of this domain seems to have no observable effect on the radiating 
shockwave within flow field. 
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Figure 49. Case 25 (Base), Time Step (a) 10 (b) 30 (c) 50 (d) 111msec 

Next, case 26 involved setting the LB conditions inside the maximum BC radius 
to observe its effect on UNDEX during the identical time steps of case 25. No 
differences were observed at time step 10msec where initial shockwave intercepts the 
water-air interface, shown in Figure 50(a). However, subsequent time steps clearly show 
distinction on the effect of wall boundary condition set inside the BC radius. First, 
Figure 50(b) shows BC at a maximum radius similar to that of the base case, but with 
shallower maximum depth. Second, Figure 50(c) shows the presence of secondary 
shockwave that is reflected off the LB as it travels back towards the charge location, 
creating turbulence and additional cavitation following its path as it intercepts the 
hammer pressure shockwaves. At this point, the additional cavitation following this 
wave is also collapsing, creating another shockwave, although smaller in magnitude. As 
this reflecting wave intercepts the left boundary condition and again reflected back into 
the water domain, the BB waves are seen disrupting this process as it begins to dissipate 
completely post time step of 111msec, shown in Figure 50(d). 
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Figure 50. Case 26 (Inside BC), Time Step (a) 10 (b) 30 (c) 50 (d) 111msec 


Case 27 shows the LB placement at the edge of the BC radius (Figure 51). The 
flow field at time step 10 and 30msec are identical to base case in both shockwave 
propagation and BC’s max radius and depth. At time step 50msec of Figure 51(c), the 
reflected LB SW begins to emerge with trails of additional cavitation forming. This LB 
shockwave is delayed due to the boundary’s placement and the magnitude seems to be 
similar to that of case 26. Lastly, Figure 51(d) shows the several BB and LB shockwaves 
intercepting to initiate additional turbulence and cavitation fonnations. 
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Figure 51. Case 27 (Edge of BC), Time Step (a) 10 (b) 30 (c) 50 (d) 111msec 

The results of LB placement well outside of the BC radius (case 28) are shown in 
Figure 52. The flow field at time step 10 and 35msec are identical to the base case in 
both shockwave propagation and max radius and depth of BC. At time step 50msec, 
shown in Figure 52(c), the delayed LB reflected wave have yet to intercept the hammer 
pressure wave generated as initial BC collapses. Once again, following closely to the 
laterally reflected wave, an initiation of additional cavitation is clearly visible in Figure 
52(c). Lastly, intercepting laterally reflected and the BB SWs are shown in Figure 52(d). 
As expected this interception takes place at a time step longer than any other cases 
discussed due to the location of LB conditions and compared to previous cases the 
additionally induced cavitation seems to have dissipated completely by this time step of 
111msec. To further analyze various UNDEX parameters numerically, a point target 
location was chosen at 100ft lateral distance and surface depth from the charge location. 
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Figure 52. Case 28 (Outside of BC), Time Step (a) 10 (b) 30 (c) 50 (d) 111msec 


2. Target Pressure Analysis 

The resulting shockwaves at the target location are analyzed for various LB 
conditions (figures 53 to 55). The first distinct spikes at approximately 25msec are initial 
compressive shocks as it reaches the target, immediately followed by tensile SW as 
water-air interface expands and spallates to initiate BC formation, the BC closure or 
“hammer” SW and followed by laterally returning reflected shockwave’s interception of 
target location. The BC collapsing shockwave is considered to be secondary and laterally 
reflected wave is considered tertiary during this case study. The other shockwaves 
observed post 110msec is BB SWs. The location and magnitude of these SWs are 
consistent theoretically. However, the timing and where both LB and BB reflected SW is 
located during the entire pressure time step can cause undue damage to the ship by 
invoking modal response. Another characteristic to point is the fact that the magnitude of 
the LB and BB shockwaves are at time even greater than that of BC collapsing SWs. 
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Target Pressure for Various Lateral Boundary Conditions 



Figure 53. Target Pressure for Various Lateral Boundaries 


Target Pressure for Various Lateral Boundary Conditions 



Figure 54. Initial, Secondary and Tertiary Target Pressure for Various Lateral 

Boundaries 
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The initial shockwave magnitude for all LB cases are identical in magnitude at 
approximately 55psi and shaped similarly; however, what follows after the BC collapse, 
differs depending on the size of the water domain or location of LB conditions. Between 
approximately 25 and 35msec where BC collapses, the magnitude of cyclic trapped 
shockwave between the surface and target location is largest for LB condition located 
outside BC max radius (Figure 54). Next, the magnitude of BC collapse shockwave also 
differs with this boundary condition at 22 psi. What follow next are three tertiary 
pressure peaks showing laterally reflected waves for all three boundary conditions with 
closest boundary showing highest and earliest magnitude of reflected shockwave. Lastly, 
the group of tensile and compressive shockwave after 110msec show various shockwaves 
as both BB and LB reflected shockwave intersect to cause additional cavitations and 
shockwaves where it reaches peak magnitudes of greater than hammer shockwave 
(Figure 55). 


Target Pressure for Various Lateral Boundary Conditions 



Figure 55. Bottom Bounce Target Pressure for Various Boundaries 


62 































3. Vertical Take-Off Velocity Analysis 

The VTO for various LB conditions are shown in Figure 56. The initial upward 
velocity peaks at 1.75ft/s for all LB conditions. What follows next mirror the pressure 
profile seen in previously. As the initial peak VTO settles to near zero velocity, the 
hammer shockwave at approximately at 37msec induces another VTO back to .25ft/s. As 
VTO reattempts to settle, the three additional peak pressure waves induced by LB 
between 49 to 75msec causes three additional peaks of VTO. The VTO settles once 
again, and then picks up approximately 120msec when collision of LB and BB reflected 
SWs cause additional BC and dynamic motion of water within the domain. It is 
important to note that VTO does not return below zero from approximately 35 to 
130msec and continues to be induced with additional forces and its velocity caused by 
BC, collisions of LB and various reflected wave. 


Vertical Take-Off Velocity (VTO) for Various Lateral Boundary Conditions 



Figure 56. VTO for Various Lateral Boundaries 


63 


























4 . Bulk Cavitation Analysis 

The BC volumes at various LB conditions compared to the base condition are 
depicted in Figure 57. Observing the initial BC formation and dissipation from time 
steps 0 to 50msec, there are several distinct differences in various LB conditions. First, 
base and LB2.3 conditions are practically identical in nature except there is a distinct 
“hump” (1.25x10 ft ) initiated at 50msec and also post 135msec. This phenomenon, also 
observed in Figure 52(c), is the initiation of additional cavitation at LB at time of initial 
shockwave reflection. This “hump” in BC volume is not observed in base case obviously 
due to different LB condition setting as “free” condition. Second, looking at the LB2.1 
and LB2.2, a similar hump is visible at an earlier time step (near 37msec) due to closer in 
proximity of the wall. Lastly, the overall BC volume of LB2.1 is much greater than rest 
compared to the base and other LB cases. Like shallower depth BB shockwave inducing 
enonnous volume of additional cavitation, the same effects are observed here for LB 
conditions well within the BC’s max radius. 


x io 5 Bulk Caviation Volume for Various Lateral Boundary Conditions 



Figure 57. Bulk Cavitation Volume for Various Lateral Boundaries 
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The next range of time step from 50 to 100msec shows the effect of LB’s 
interferences and additional cavitation formed as the lateral shockwave intercepts the 
hammer shockwave. Clearly, the magnitude and speed of the additional cavitation 
formed are much higher for LB2.1. The post 100msec ranges show drastic increase in 
additional cavitation fonnation due to BB inclusion of already chaotic water domain. 
Consistent in behavior to the 50 to 100msec range, the LB that is the closest to the charge 
had the largest additional cavitation formed surpassing the initial cavitation volume for 
all cases. It is surprising to note that the LB conditions have almost the same effect in 
inducing additional cavitation as the BB has on BC. 

B. VARYING REFLECTIVITY INSIDE THE BC RADIUS 

In order to further investigate the degree of shockwave propagation as it pertains 
to the percentage of reflectivity of the LB conditions, case 26 (LB2.1) were simulated 
once again using 50% and free boundary conditions. For case 27 (LB2.2) and 28 (LB2.3) 
the effects of such change in reflectivity due to greater distance, in reference to max BC 
radius, were not as drastic as case 26. All other parameters were maintained as 
previously described in Table 3. 

1. Flow Field Analysis 

The flow field results for 50% reflective LBs are shown in Figure 58. 
Comparable to base case 26, where boundary condition was set as 100% reflective (wall), 
the general features of radiating initial shockwave, BC fonnation, reflected shockwave 
and its interception of BB waves and BC collapse wave are nearly identical. Case 26 is 
the base case in this section where its LB conditions is set as wall and case 26.1 is set at 
50% reflectivity. As expected, the magnitude of shockwave for the 50% reflective wall is 
clearly lower than the case 26 but its shockwave travels at a pattern consistent with case 
26. The flow field at 111msec is drastically different than that of base case 26. As 
Figure 58(d) shows, the magnitude in which BB wave interacts with laterally reflected 
waves are a lot higher than the base case. Since, 50% of laterally reflected energy was 
absorbed by the LB; the intercepting wave has less energy to neutralize BB waves as it 
traveling up to the water-air interface. In short, as BB waves travels up and passes the 
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charge depth, the energy dissipation is a lot lower for 50% lateral wall reflective case 
26.1. The result is increased UNDEX forces applied to object in its path and on surface 
of water. 



Figure 58. Case 26.1 (Inside of BC, 50% Reflectivity), Time Steps (a) 10 (b) 30 

(c) 50 (d) 111msec 

Next, the flow field result for Case 26.2 is shown in Figure 59. Case 26.2 depicts 
case 26 with free or no reflective LB conditions well within the max radius of BC. The 
flow field result is almost identical to the base case 25 for all chosen time step, except for 
111msec. At approximately 100ft deep below charge depth there is a pressure spike 
where BB SWs meet the downward traveling BC hammer SW (Figure 59(d)). As seen, 
where these two pressure waves collide, there is a noticeable spike that is similar to when 
the laterally reflected wave and BB wave intercept. Although this seems minor, the 
distinct and unexpected difference in pressure spike can cause additional forces applied to 
the floating object within the timeframe of UNDEX event. 
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Figure 59. Case 26.2 (Inside of BC, No Reflectivity), Time Steps (a) 10 (b) 30 


(c) 50 (d) 111msec 

2. Target Pressure Analysis 

The pressure profde for cases 25, 26, 26.1 and 26.2 are shown in figures 60 to 62. 
Observing these pressure profiles of four different cases, the tertiary LB SW following 
BC hammer shockwave is clearly absent except for case 26 and 26.1 where LBs were set 
as a rigid wall and 50% reflective wall. The observed shockwaves are expected for both 
cases in terms of their location and magnitude. A closer look at timeframe from 20 to 
50msec (Figure 61) shows base case 25 and case 26.2 pressure profiles are indeed 
identical for both BC collapse shockwaves of 17.5psi and absence of LB reflection spike 
near 48msec. Likewise for case 26 and 26.1, the pressure profiles also closely resemble 
each other, being identical in BC collapse pressure spike of 22psi at 37msec and presence 
of laterally reflected wave near 48mse, with wall LB setting having the higher pressure 
spike at 23psi. 
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Target Pressure for Case 26 at Various Reflectivity 



Figure 60. Target Pressure for Varying Boundary Reflectivity 


Target Pressure for Case 26 at Various Relectivity 



Figure 61. Initial and Secondary Target Pressure for Varying Boundary 

Reflectivity 
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Taking a closer look at time step ranges from 100 to 150msec where BB SWs 
start to emerge, again the presence or lack of LB conditions play a crucial role in 
invoking additional shockwaves during this time frame. Both cases 25 and 26.2 show 
similar pressure fluctuations between 13.5 and 17psi. Cases 26 and 26.1 show pressure 
fluctuations between 6 and 24psi. The differences between wall and free LB conditions 
are distinct with possibility of LB condition invoking higher pressure fluctuations and 
within the changes, the reflectivity setting closely resembles how much fluctuations are 
observed for cases 26 and 26.1. 


Target Pressure for Case 26 at Various Reflectivity 



Figure 62. Bottom Bounce Target Pressure for Varying Boundary Reflectivity 

3. Vertical Take-Off Velocity Analysis 

The VTO results are shown in Figures 63. The initial spike in VTO peaks at 
1.75ft/s for all cases. The secondary or LB induced VTO for all cases are different to 
some degree, but the presence of spike is apparent for presence of LB setting at 50% or 
wall. For cases 25 and 26.2, the spike is absent; however, there is a minor difference 
between the two. Again for BB induced VTO, there are distinct differences between 

69 































presence of LB conditions and higher magnitudes are observed for cases where 50% or 
wall setting is present. Lastly, compared to cases 27 and 28 (Figure 56), the absence of 
varying LB induced spikes, as well as, reduced fluctuation in BB induced VTO are also 
apparent for this case (Figure 63). 


Vertical Take-Off Velocity (VTO) for Case 26 at Various Reflectivity 



Figure 63. VTO for Various Lateral Boundary Reflectivity 

4. Bulk Cavitation Analysis 

Taking flow field results into considerations, the BC volume formations were also 
analyzed numerically and shown in Figure 64. The changing BC volume of case 26.1 
and 26.2 were compared to base case 25 and case 26 (Figure 57). There are two 
distinctly unexpected characteristic observed for cases 25 and 26.2 (free) vs. 26 and 26.1 
(50% or wall). As mentioned, case 25 is base case with free boundary and case 26.1 is 
the same except the location of the LB. In case 25, the boundary is located well outside 
the BC max radius and case 26.1 the boundary is located inside the BC max radius. The 
volumetric formation of BC result shows that results of case 25 and 26.2 are not identical. 
During the initial collapse of BC near 40 to 50msec, the results show case 26.2 to 
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collapse a lot faster than that of case 25. Additionally, near 30msec and beyond, the case 
26.2 shows greater volume spike than case 25. This is during the secondary pressure 
spike where BC hammer is present (approximately 35msec) and from the observed 
results, case 26.1 and 26.2 had the most impact in their BCs. For time step ranges of 50 
to 100msec, where LB reflective shockwaves are present, case 26 and 26.1 had the most 
impact. This was expected due to presences of reflective setting in their LB conditions. 
Lastly, time step ranges from 125 to 150msec show cases 26 and 26.1 with most 
additional cavitation formed due to same reasoning behind presence of reflective LB 
conditions. 



Figure 64. Bulk Cavitation Volume for Various Lateral Boundary Reflectivity 
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VII. PRECURSOR TO FULLY COUPLED RUNS 


A. THE BLOCKED CELL METHOD 

The presence of foreign body and its effect on UNDEX parameters in LOD is of 
interest for this chapter. Likewise, UNDEX parameter’s impact on the FSA is also of 
interest. As a precursor to including rigid body or Lagrangian solid in the Eulerian 
domain for fully coupled run, use of “blocked cell” method was utilized. The blocked 
cell method is set and invoked during Pregemini phase, where inclusion of such object in 
the Eulerian domain acts as a rigid body in lieu of an actual Lagrange structure in its 
place, without actually running a fully coupled run. The placements of foreign bodies in 
three locations as blocked cells (BK) are shown in Figure 65. Generality exists to show 
that placements of foreign surface object near origin, at midpoint and at the edge of BC 
radius will have an observable impact on UNDEX parameters; hence, making these cases 
worth the look prior to fully coupled simulations. As before, the target, charge and 
DYSMAS set up are as described in Table 4. 


A Y 



Air Column 


CHARGE 

(X,Y,Z) C 


Water Column 



Figure 65. Blocked Cell Methods 
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Case 

# 

Run 

Mode 

Euler 

Dimesions 
(X, Y, Z) [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

29 

2D, Euler 

300, 1, 300 

10 

15, 105, 255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

30 

2D, Euler 

300, 1, 300 

10 

15 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

31 

2D, Euler 

300, 1, 300 

10 

105 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

32 

2D, Euler 

300, 1, 300 

10 

255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 


Table 4. Case Studies for Blocked Cell Methods 


B. SUMMARY OF RESULTS AND ANALYSIS 

1. Llow Lield Analysis 

The flow field results for cases 29 through 32 are depicted in figure 66 through 
69. The base case, without any blocked cell is shown in Figure 66. No unusual 
observation was made and all time steps with significant UNDEX events following the 
same pattern as previous cases studied: (a) initial SW propagation (b) BC reaching 
maximum radius and depth (c) collapse of BC and initiation of hammer SWs (d) 
emergence of BB SW as it travels to intercept water-air interface to invoke additional 
cavitation. Two things that stand out are the characteristics of residual cavitations shown 
in Figure 66(d). 



Figure 66. Case 29 (Base), Time Step (a) 10 (b) 30 (c) 50 (d) 100msec 
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Case 30 with blocked cell closest to the detonation is depicted in Figure 67. The 
distinctly different characteristics are (a) observation of local cavitation near and around 
the blocked cell, (b) disruption in BC formation by the blocked cell, (c) presence of 
additional hammer SWs originating from collapse of local cavitation of blocked cell and 
speed at which initial BC collapses. 



Figure 67. Case 30 (BK1), Time Step (a) 10 (b) 30 (c) 50 (d) 100msec 


Case 31 with blocked cell at midpoint to BC radius is depicted in Figure 68. The 
distinctly different characteristics are (a) lack of local cavitation near and around the 
blocked cell, (b) disruption in BC formation by the blocked cell as it fully matures and 
local cavitation, (c) presence of additional hammer SWs originating from collapse of 
local cavitation of blocked cell merging with the initial BC collapse and (d) emergence of 
laterally reflected SW originated from the local cavitation hammer shockwave. 
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Figure 68. Case 31 (BK2), Time Step (a) 10 (b) 30 (c) 50 (d) 100msec 


Case 32 with blocked cell located well outside the max radius of BC is depicted in 
Figure 69. Once again, the distinctly different characteristics are (a) lack of local 
cavitation near and around the blocked cell, (b) lack disruption in BC formation and 
presence of local cavitation, (c) lack of additional hammer SWs originating from collapse 
of local cavitation and (d) emergence of laterally reflected SW originated from unknown 
source. The origin of this laterally reflected wave can be speculated from local cavitation 
collapse from blocked cell, but with time step and location of blocked cell for the case 32, 
this is unlikely. 
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Figure 69. Case 32 (BK3), Time Step (a) 10 (b) 30 (c) 50 (d) 100msec 
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2. Target Pressure Analysis 

The target pressure profiles for each blocked cell location are depicted in figures 
70 through 72. Figure 70 is target pressure profile for BK1 compared to the base case 
without any blocked cells. As can be seen, although there are no impacts to the 
magnitude of initial shockwave of nearly 700psi near 10msec, a significant differences in 
BC collapse shockwave’s time and magnitude can be observed. While the base case’s 
BC collapsing shockwave reaches close to 400psi, the BK1 only experiences close to 
lOOpsi. BK1 also experiences BC closure almost 10msec faster than base case. As 
observed in the flow field results, the presence of BK1 also expedites the BC collapse. 


Pressure Profile vs. Time 



Figure 70. Pressure Profile for BK1 and Base Case 


Target pressure profile for BK2 and base case is depicted in Figure 71. Once 
again, the initial shockwaves are almost identical in magnitude and timing with close to 
240psi at time 25msec. Unlike previous case, this time the BC collapse shockwave is 
higher for the BK2 at close to 85psi while base case hovers at 25psi. Post 100msec show 
presence of additional cavitations formed when shockwaves originating from the collapse 
of local cavitation is reflected off the lateral boundaries, shown in flow field of Figure 68. 
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Pressure Profile vs. Time 



Figure 71. Pressure Profile for BK2 and Base Case 


Pressure Profile vs. Tme 



Figure 72. Pressure Profile for BK3 and Base Case 


Target pressure profile for BK3 and base case is depicted in Figure 72. Once 
again, the initial shockwaves are almost identical in magnitude and timing with close to 
lOOpsi at time 46msec. Similar to previous case, this time the BC collapse shockwave is 

higher for the BK3 at close to 35psi while base case hovers at 20psi. Post 120msec also 
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show presence of additional cavitations fonned when shockwaves originating from the 
collapse of local cavitation is reflected off the lateral boundaries, shown in flow field of 
Figure 69. 


3. Vertical Take-Off Velocity Analysis 

The VTO for each blocked cell location are depicted in figures 73 through 75. 
The VTO for base case and BK1 are depicted in Figure 73. Both cases’ VTO peaks at 
11 ft/s followed by an immediate collapse and fluctuation near and around 3ft/s. The 
VTO fluctuations for the base case seems to be more dramatic in that lack of blocked cell 
that can hinder movement will freely let water to slush around the domain. Conversely, 
the presence of blocked cells greatly decreases the movement of the water around it since 
the energy is transferred to and from blocked cell invoking a beginning of FSA 
phenomenon. 
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Figure 73. VTO for BK1 and Base Case 

The VTO for the base case and BK2 are depicted in Figure 74. Both cases’ VTO 

peaks at 1.5ft/s followed by an immediate collapse and fluctuation near and around 

0.25ft/s. The VTO fluctuations for the BK2 seem to be more dramatic with dip of -1.5ft/s 
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Vertical Take-Off Velocity (VTO) vs. Time 
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at approximately 40msec. The dip at 40msec is counterintuitive since the blocked cell is 
not moving, but further investigating flow field results (Figure 68), show that this is due 
to traveling shockwave originated from the collapsing local cavitation. 


Vertical Take-Off Velocity (VTO) \s. Time 



Figure 74. VTO for BK2 and Base Case 


Lastly, VTO for base case and BK3 are compared in Figure 75. For this case, the 
presence of initial shockwaves at 0.3psi at 55msec and presence of shockwave 
originating from collapse of local cavitation are close to being identical, except, for BK3, 
the VTO continues to increase throughout the time step. Since blocked cells are treated 
as a rigid material and reflecting surface, BK3’s characteristic can only be explained as 
error in Gemini’s solution [12]. Numerous attempts were made to fix this error through 
remeshing Eulerian domain and placement of the blocked cells, as well as, target 
position. At this time, an acceptable solution does not exist although the pattern and 
presence of initial and additional shockwave are consistent with previous cases. 
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Vertical Take-Off Velocity (VTO) vs. Time 



Figure 75. VTO for BK3 and Base Case 


4. Bulk Cavitation Analysis 

The BC results for all blocked cell cases are depicted in Figure 76. Similar 
characteristics to previous cases are still observed, where BC is initiated at 10msec when 
initial shockwave reaches water-air interface, reaches its maximum radius at 
approximately 35msec, dissipates completely by 60msec and additional cavitation 
observed near 120msec where BB reaches water-air interface to invoke this phenomenon. 
Two distinctly different patterns emerge for BK2 where the BC starts to collapse 22msec 
with another increase near 40msec when local cavitations are invoked and completes 
dissipation near 60msec. For BK3, lack of local cavitation reduces the overall BC 
volume; hence, its collapse prior to 60msec can also be observed. All four cases 
experienced additional cavitation formations starting at 120msec when BB shockwave 
emerged to initiate another set of BC near water-air interface. 
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Bulk Caviation Volume vs. Time 



Figure 76. Bulk Cavitation Volume for Blocked Cell Cases 


82 



























VIII. CONCLUSION AND RECOMMENDATION FOR FURTHER 

STUDIES 


A. CONCLUSIONS 

An in-depth analysis and characterization of littoral ocean domain (LOD) were 
conducted for transient phase of UNDEX by analyzing resulting parameters such as 
initial shock wave, vertical take-off velocity and bulk cavitation. Implementing various 
conditions such as ocean depth, charge size and depth, boundaries for bottom and side 
conditions of fluid domain and even including a small rigid body, these results were 
compared and contrasted to the current Eulerian fluid model. 

First, results of varying ocean depths show clear distinction in UNDEX 
characterization at depth shallower than 300ft. At this depth, multiple shockwaves and 
resulting vertical take-off velocity exists due to additionally induced cavitation resulting 
from the reflecting pressure wave from the shallower bottom. The reflecting shockwave 
and energy contained within are the main culprit of shallow UNDEX characteristics. 

Second, the results of varying charge size and depth showed that charge size of 
less than or equal to 3001bs showed linear relationships and as the charge depth reached 
water-air or water-bottom interface, the characteristics of resulting UNDEX parameters 
became increasingly amplified and chaotic. The combined bottom bounce and initial 
shockwave for charge depth at water-bottom interface showed combined effects of its 
resulting UNDEX parameters with greater BC volume while delayed initial response. 

Third, results of varying lateral boundary conditions show that as the lateral 
boundary distance is brought closer to the detonation source, inside the radius of bulk 
cavitation, its behavior and the rest of the UNDEX parameters also become increasingly 
amplified due to similar effects seen in shallower bottom depths. Varying its reflectivity 
also amplified UNDEX reactions; however, differences in reflectivity between free and 
50% lateral boundary condition inside and outside the radius of BC showed negligible 
changes in the UNDEX results. 
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Lastly, adding blocked cells prior to a full scale coupled run showed that fluid 
behaves more erratically as these blocked cells are situated within the radius of forming 
BC. As fluid behaves differently around the existing small rigid blocks, so too are the 
UNDEX parameters due to reflective energy from these rigid blocks, invoking energy 
back into the fluid domain. As a result, a precursor to how fluid may behave starts to 
emerge for a fluid structure interaction or fully coupled run cases. 

As a result, recommended characterization and boundaries of LOD for UNDEX 
studies start to emerge. From this study, the ocean depth of 300ft, lateral boundary 
distance of greater than 1/3 radius distance from the edge of bulk cavitation, charge 
weight of 3001bs or less and its placement at minimum of 1/3 distance away from the 
existing boundaries, such as bottom or lateral wall are recommended. For fully coupled 
runs, the placement of Lagrangian solid within the bulk cavitation radius and inclusion of 
various bottom sediments to increase realism with caveats of setting its reflectivity 
condition to rigid or wall are also recommend for LOD. Although cyclic bubble was not 
studied extensively during this research, the obvious choice of bubble setting is to pick a 
charge size that will produce an initial bubble size consistent with previous LOD 
boundaries. If purpose is to study the transient phase of UNDEX then this is not a huge 
factor. However, if whipping analysis is of interest, choosing a charge size that will 
invoke a bubble radius of at minimum of less than chosen LOD depth is recommended. 

B. RECOMMENDATION FOR FURTHER STUDIES 

While this study shed some light on the benefits of establishing the boundaries of 
LOD for future coupled runs of various Lagrangian models or an actual LCS FE model, 
due to the sheer size of the problem in computing time and efficiency, other parts of 
Eulerian characteristics such as implementing viscous code, adding shelves within the 
fluid domain to study the effects of its angle and further study of steady state UNDEX 
phenomena like cyclic bubbles and its effect are recommended. 
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APPENDIX A. COMBINED DYSMAS/GEMINI SIMULATED CASES 


Case 

# 

Run 

Mode 

Euler 

Dlmesions 

(X,Y,Z) [ft] 

Target 

Depth 

[ft] 

Target 

Distance 

[ft] 

Charge 

Type 

[ft] 

Charge 

Weight 

[lbs] 

Charge 

Depth 

[ft] 

Bottom 

Type 

B1 

[zmax] 

B2 

[xmax] 

B3 

[zmin] 

B4 

[xmin] 

Grid 

Size 

[ft] 

Remarks 

1 

2D, Euler 

800, 1 , 1000 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying Ocean Depth from 1000ft to 75ft 
deep. All other UNDEX parameters, 
including charge size and charge depth 

maintained constant. 

2 

2D, Euler 

800, 1, 900 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

3 

2D, Euler 

800, 1, 800 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

4 

2D, Euler 

800, 1, 700 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

5 

2D, Euler 

800, 1, 600 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

6 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

7 

2D, Euler 

800, 1 , 400 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

8 

2D, Euler 

800, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

9 

2D, Euler 

800, 1, 200 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

10 

2D, Euler 

800, 1, 175 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

11 

2D, Euler 

800, 1, 150 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

12 

2D, Euler 

800, 1 , 125 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

13 

2D, Euler 

800, 1 , 100 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

14 

2D, Euler 

800, 1, 75 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

15 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying Charge Size. All other UNDEX 
parameters maintained constant. 

16 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

200 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

17 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

300 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

IS 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

400 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

19 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

500 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

20 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

100 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Varying Charge Depth. All other UNDEX 
parameters maintained constant. 

21 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

200 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

22 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

300 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

23 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

400 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

24 

2D, Euler 

800, 1, 500 

0.820 

100 

HBX-1 

100 

500 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

25 

2D, Euler 

230, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

Base Case, free boundary @ xmax 

26 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

Varying Case, wall boundary @ xmax. For 
26.1 & 26.2, varied reflectivity 50% &free 
for boundary located inside BC zone. 

26.1 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

0.50 

Wall 

Wall 

0.656 

26.2 

2D, Euler 

164, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.656 

27 

2D, Euler 

197, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

28 

2D, Euler 

230, 1, 300 

0.820 

100 

HBX-1 

25 

50 

Clay 

Free 

Wall 

Wall 

Wall 

0.656 

29 

2D, Euler 

300, 1, 300 

10 

15, 105, 255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Base Case, No Blocked Cell 

30 

2D, Euler 

300, 1, 300 

10 

15 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked cell at position (1) 

31 

2D, Euler 

300, 1, 300 

10 

105 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked cell at position (2) 

32 

2D, Euler 

300, 1, 300 

10 

255 

HBX-1 

25 

50 

Clay 

Free 

Free 

Wall 

Wall 

0.164 

Blocked cell at position (3) 
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APPENDIX B: BULK CAVITATION THEORY 


The following section describes how the method of Arons is used to map out the 
bounds of the cavitated region that forms when surface cutoff tries to lower the absolute 
pressure below the cavitation pressure. A continuation from section II. B. 2 of this study, 
it is a direct summary from the critical sections within the technical paper published by 
Costanzo and Gordon in May 1983, titled “A Solution to the Axisymmetric Bulk 
Cavitation Problem” and for that I am grateful for its use and many thanks to Contanzo 
and Gordon for use of their work as primary reference in characterizing BC for this study 

[9]. 

A. DERIVATION OF METHODS OF ARONS AND TANGENT RULE 

The target pressure (P) at its location of (X, Y) is assumed to be the absolute 
pressure prior to arrival of the reflected wave from water-air interface. This pressure is 
culmination of overpressure generated by SW, atmospheric and hydrostatic pressure. 
Now, let P be function of radial (r) and angular (a) coordinate originating from the image 
charge, W;, as shown in Figure 20. The upper cavitation boundary is defined as the locus 
of points at which the absolute pressure drops to the cavitation pressure upon arrival of 
the rarefaction wave [9]. At this point, water at cavitation will remain cavitated as long 
as the absolute pressure does not rise above the vapor pressure of water. For all practical 
purpose for the method of Arons, the vapor and cavitation pressure are assumed to be 
zero. Thus, the equation defining the upper cavitation boundary is 


P[a,r^j + P r = 0 


(B.l) 


If P r , the reflected wave, is expressed in terms of the charge weight and standoff 

1 /3 B 

using the method of images, then P r = -A(W /r) . Substituting this expression into 
Equation (B.l) gives 


P(a,r)-A 




= 0 


v r j 
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(B.2) 



Where, A and B are constants specific to the charge or explosive types. Now, before the 
cavitation boundary is discussed, the tenn breaking pressure must be defined. The 
breaking pressure is the magnitude of the rarefaction wave (or relief wave, as commonly 
tenned) which reduces the absolute pressure at a point to the cavitation pressure. In other 
words, since the cavitation pressure is taken to be 0 psi absolute, the breaking pressure 
has a magnitude equal to the absolute pressure at a point prior to the occurrence of 
cavitation at that point. The equation of the lower cavitation boundary is derived from 
consideration of the propagation of this breaking pressure into the uncavitated water 
beneath this boundary. Let P r = P(a, r) be the breaking pressure for a point lying at the 
lower cavitation boundary. At this lower boundary, P(a, r) must propagate into 
uncavitated water with spherical attenuation resulting in P r = P(a, r)(R/r) B . This is 
represented in Figure 77. 



Figure 77. Propagation of Breaking Pressure into Uncavitated Water, from [9] 

This pressure expression must have a faster decay rate than the general expression 
for absolute pressure, P = P(a, r), or the water located along the dashed line extending 
below the lower cavitation boundary will continue to cavitate. Thus, at a point on the 
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lower cavitation boundary, the decay rate of both of these pressure expressions is the 
same. Therefore: 


±[ P ( a ,r)(Rlryl^[P(a,r)l 


dr 1 


(B.2.a) 


If all terms of this equation are shifted to one side of the equal sign and r’ B is 
factored out, then application of the chain rule to differentiation of a product results in: 


d_ 

dr 


[r B )p(a,r)-P(a,r){R B ) 


r=R 


0 


(B.2.b) 


This expression may be simplified further by recognizing that the product, 
P(a,r)(R B ), pertains to a specific point and thus may be treated as a constant when 
differentiating. Therefore, this equation becomes: 


= 0 

(B-3) 

where r = R, the point at the lower cavitation boundary for a given a. Equation (B.3) is 
the method of Arons for determining the lower cavitation boundary. 

In Figure 78, the general shape of the BC envelope as determined by the method 
of Aron is given. It is of fundamental interest to the BC problem to locate the point at 
which a line drawn from the image charge is tangent to the upper cavitation boundary. 
For the derivation of this tangent point, both sides of Equation (B.2) are multiplied by r B 
and the resulting equation and its total derivation with respect to r are written below as 
equations (B.4) and (B.5), respectively. 


— (r‘)p(a,r) 


[r B )p(a,r)-A[w m ) B =0 


(B.4) 


dr L 


[r B )p(a,r)-A[w^) L 


' _ d 

\(r B )p(a,r)\ 

d 

+ — 

r(r B )R(a,r)] 

r da^ 

~~dr 

LV / V 'A 

da 

LV / v 'A 

ydr , 


= 0 


(B.5) 
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For a given value of da/dr, the simultaneous solution of these two equations gives 
the corresponding values of a and r. When the ray extending from the image charge is 
tangent to the upper cavitation boundary, a is a maximum and da/dr = 0. When this 
value for da/dr is substituted into Equation (B.5), the following is obtained: 


j-[(r’)P( a ,r) 


0 


(B-6) 


Therefore, the simultaneous solutions of Equations (B.4) and (B.6) give the value 
of a and r at the tangent point. Equation (B.6) also happens to be equivalent to Equation 
(B.3), the equation whose solution determines the lower cavitation boundary. This 
indicates that the lower and upper cavitation boundaries intersect at the point at which a 
line drawn from the image charge is tangent to the upper boundary. This rule of tangency 
is also illustrated in Figure 78. 



Figure 78. Bulk Cavitation Bounds and Rule of Tangency, from [10] 

B. DEVELOPMENT OF THE CLOSURE MODEL 
1. General Description 

The general representation of a point which lies within the cavitated region is 
given in Figure 79. Upon the arrival of the relief wave at this point, the vertical 


component of the instantaneous water particle velocity is the vector sum of the vertical 
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components of the two velocity vectors as shown. This vector sum is termed the vertical 
water particle kickoff velocity or vertical take-off velocity (VTO) at point (X, Y) and is 
dependent solely upon the magnitude and directions of the incident and reflected acoustic 
water particle velocities. The directions of these velocities are defined by the unit 

vectors, 1 and J , as shown in the figure. 



At every horizontal range within the extent of the BC envelope, there exists a 
column of water which lies between the surface and the upper cavitation boundary. Since 
this water does not cavitate, it is assumed that all water particles contained in this vertical 
column are kicked off simultaneously with kickoff velocity of the water particle located 
at the upper cavitation boundary. This is illustrated in Figure 80. This surface layer, 
which is kicked off at the time of relief wave arrival at the upper cavitation boundary 

with an initial velocity, , is decelerated by atmospheric pressure and gravity. At a short 
time later, the water particle lying directly beneath this surface layer is kicked off with an 

initial velocity, 0 . 

Since this particle lies within the cavitated region and is kicked off after the 

surface layer, it becomes separated from the surface layer and therefore has no 

atmospheric pressure acting upon it. Thus, only gravity decelerates this particle. 

Eventually, this water particle will collide with the surface layer above it. In the 

91 






development of the closure model, it is assumed that this is a perfectly inelastic collision. 
Thus, this particle and the surface layer above it now form an augmented surface layer 
which has a velocity derived from the conservation of momentum consideration. As with 
the original surface layer, this augmented surface layer has atmospheric pressure and 
gravity acting to decelerate it. 



Figure 80. Kickoff Velocity of the Surface Layer, from [9] 

Since the particles lying below the original surface layer are all kicked off at 
different times with different vertical kickoff velocities, inelastic collisions will occur one 
at a time between the growing surface layer and the particles directly below it. This 
process is known as accretion. If the surface layer displacement history at a horizontal 
range, X, is plotted with t=0 referring to the time of explosives charge detonation, the 
curve in Figure 81 is obtained. This curve is not quite a perfect parabola due to the fact 
that the surface layer mass is changing. Also note that this curve accounts for the initial 
displacement of the surface layer due to the incident SW prior to relief wave arrival at the 
upper cavitation boundary. 
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Figure 81. Surface layer Displacement, from [9] 


The vertical component of the water particle velocity for the point which lies at 
the lower cavitation boundary at this same horizontal range, X, also must be detennined. 
Since cavitation does not extend below this point, separation will not occur between the 
underlying water particles. Hence, the water which lies below this point at the lower 
cavitation boundary will have a vertical velocity component which is dependent upon the 

varying velocities ^ and ^ and their corresponding afterflow terms. This is depicted for 
point (X, Y l ) in Figure 82. 



Figure 82. Velocity of a Point at the Lower Cavitation Boundary, from [9] 
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2. Bulk Cavitation Closure Pulse 

The magnitude of the pressure pulse generated by the collision of the accreting 
lower cavitation boundary at closure will now be detennined. Since in the closure model, 
the vertical components of the velocities of both the surface layer and the lower 
cavitation boundary for a particular horizontal range, X, are known at the time of closure, 
the closure pressure pulse magnitude, P c , can be calculated by multiplying one half the 
relative velocity of impact of these two water masses by the characteristic impedance of 


the medium. That is: 


Pc = ^(V L -V u ) 


(B.7) 


where V is positive upward and Vu and Vl refer to the vertical velocity components at 
closure of the surface or upper layer and the lower cavitation boundary, respectively. The 
closure pressure pulse consists of two compressive waves of magnitude, Pc; one that 
travels upward and one that travels downwards. These closure pressures are most readily 
superimposed upon the incident SW and reflected pressure at the horizontal range at 
which closure first occurs, since at the range, the direction of propagation of the closure 
pulse is vertical. At other horizontal ranges, the direction of propagation of the closure 
pulse varies from the vertical direction by an angle vp, as shown in Figure 83. According 
to the work done by Cushing [10], if C is the sound speed in water and Vc is the speed at 
which the cavitated region is closing at some horizontal range, X, then: 


sin y/ = 


(B.8) 


Thus, the magnitude of the pressure pulse which propagates along a path oriented 
at this angle vp is determined by dividing the expression of Equation (B.7) by cos ip. That 


pc<y L -v v ) 

2 cos y/ 


(B.9) 
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Figure 83. Closure Pulse at a Horizontal Range Different from that of First 

Closure, from [9] 


The duration time of the closure pulse at the depth of closure is determined by the 
method of images and is expressed by: 


-cos y/ 


(B-10) 


where Dc is the depth of closure at the horizontal range of interest, X. Thus, regardless 
of the angle of propagation of the closure pulse, the impulse (pressure X duration) 
associated with cavitation closure at a particular horizontal range is the same and may be 
expressed as shown in Equation (B.10). 


Im pulse = p(V L - V u )D C 


(B. 11) 
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APPENDIX C. BULK CAVITATION CASE STUDIES 


Following case studies were conducted in order to better characterize behaviors of 
bulk cavitation prior to its inclusion into shallow water or littoral ocean domain 
environments. The MATLAB codes used for this study is included in the Appendix of 
this report. 

1. Comparison of 2-D BC Zone, varying charge type, charge weight and charge 
detonation depth (Appendix C contains MATLAB codes for 2D BC Zones): 

a. 200 lb charges of HBX-1, TNT and PETN at a depth of 25 ft : 

Figure 84 is a 2D representation of the bulk cavitaion zone of a 2001bs charge 
denotated at a depth of 25ft. The charge lies along the left vertical edge of each subplot. 
Each subplot represents a different explosive type as titled. This figure was created for 
the sole purpose of comparing the BC zone of the three different types of explosive with 
the other variables (depth of charge, weight of charge) held constant. 



Figure 84. 2-D Bulk Cavitation of 2001bs of FIBX-1, TNT and PETN @ 25ft 

Since water cannot support a tension (i.e. negative pressure), cavitation occurs at 
points in the water column where the sum of the pressures at that point are less than or 
equal to zero. This pressure is sum of atmospheric, water depth, incident shock pressure, 
and rarefaction wave pressure. As can be clearly seen from the above figure, even 
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though all the explosive charge weights and detonation depths are the same, each 
explosive charge type creates a different BC region in both depth and radial distance 
outward. This is due to the fact that they have different constants (Ki, K 2 . A|, A 2 ) 
associated with them that are used in the calculation of their maximum pressure and their 
decay constant. 

The Pentolite (PETN) ultimately had the lowest K 2 and A 2 thus causing it to have 
the deepest but shortest BC zone. TNT and HBX-1 are very similar in their K 2 value 
however TNT’s A 2 value is significantly smaller than HBX-1, causing it to “reach” 
further outward from point of detonation. In addition, the difference in Ai and A 2 values 
seems to be the major driving factor for how far out radially the BC zone expands. The 
Pentolite had the largest difference and thus the smallest radial distance, while TNT has 
the smallest difference and thus the largest radial distance of the BC zone. 

b. 100 lb, 200 lb, and 300 lb charges of HBX-1 at 50 ft: 

Figure 85 is a 2D representation of the bulk cavitaion zone of a HBX explosive 
charge of varying weights (1001b, 200 lb, and 3001b) detonated at a constant depth of 50 
feet. Once again the location of explosive depicted in the left vertical axis of each 
subplot. The purpose of the plot is to represent the effect of charge weight on the bulk 
caviation zone with other variables (charge type, detonation depth) held constant. As 
expected, as the charge weight is increased while the detonation depth and charge type 
are held constant, both the depth and breadth of the BC zone are increased. 



Figure 85. 2-D Bulk Cavitation of 100, 200, and 3001bs of HBX-1@ 50ft 
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c. 250 lb charge of TNT detonated at varying depths ( 5, 50 and 500ft): 


Figure 86 is a 2D representation of the bulk cavitaion zone of a TNT explosive 
charge of 2501bs detonated at varying depths of 5, 50, and 500 ft. Figure 87 shows 
extended range of 250 lb at 500ft for obvious reason. Once again, the explosive is 
depicted along the left vertical axis of each subplot. Again, the purpose of the plot is to 
represent the effect of charge weight on the bulk caviation zone with other variables 
(charge type, detonation depth) held constant. 
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Figure 86. 2-D Bulk Cavitation of 2501bs TNT detonated @ 5, 50 and 500ft 
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Figure 87. Extended View of 2501bs TNT detonated @ 500ft 


As expected, as the depth of the detonation point of a charge varies, with the 
weight and type of charge being held constant, the BC region grows both in depth and 
breadth, with the exception of the 500’ depth charge. This charge does have massive 
radial distance (more than twice the radial distance of charge at 50ft) compared to the 
other two and the thickness of the BC zone seems to grow radially rather than maximum 
thickness being at or near the center of the BC zone. This is most likely due to the fact 
that the pressure from incident and image charges are very low (having decayed 
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significantly already as it traveled vertically 500ft) and are at a relatively shallow depth. 
Therefore, the total pressure at this water column, although significantly lower than other 
two depth cases, is enough to create a cavitation while the surrounding pressure rises to 
the vapor pressure of the water (about 0.3 psi). In addition, the SW that causes the 
pressure havoc in this water column and causes the BC has traveled significantly farther 
than shallower depth charges. 

2. 3-D visualization of the BC space for an underwater explosion resulting from the 
detonation of a l,0001bs charge of TNT located at a depth of 25ft: (Appendix C 
contains MATLAB codes for 3D TNT BC Zone): 

Figures 88 and 89 depict the 2D and 3D representation of l,0001bs charge of TNT 
located at 25ft depth. The growth and shape of the BC zone follows the same pattern 
observed previously in section 1 of this report. Figure 5 below shows regional and line 
depiction of the BC region. Following the same trend of BC formation, the maximum 
thickness appears at 200ft radial distance from “ground zero” or point of detonation with 
about 40ft thic kn ess. The radial length of the BC region is about 650ft. 
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Figure 88. SPY and Line View of l,0001bs TNT at 25ft Depth 

Figure 89 is a 3D representation of BC. The overall shape is almost like a flying 
saucer with a dip on the bottom side of the disc. For the below representations, the top of 
the BC region was removed to show the contour shape within the region. The figure on 
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the right depicts “ground zero” and shows the radial expansion of the explosion below the 
BC zone. The usual maximum dip that is represented at the center of the BC zone is also 
depicted in both figures almost like an upside down bell. 



Figure 89. 3D General Top and Peak Top View of lOOOlbs TNT at 25ft Depth 
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APPENDIX D. BULK CAVITATION MATLAB CODES 


Computation of Bulk Cavitation Zone 
for Underwater Explosions 
By 

LT Dean Ahn, USN 


% This program is used to develop the BC envelope for an 

% underwater charge of TNT at various 1000 lbf charge weights, 25ft depth 
% and plot them in 2D & 3D. 


%% Prep Matlab workspace 

clear all; 

clc; close all ; 

% 2-D SPY: BC visualization of UNDEX, 1000 lb charge TNT @25ft 

clear all; 

clc; 

% Constants 

Pa = 14.7; % Atmpospheric pressure (psi) 

Gamma = 62.5/144; % Weight density of water (lb/ft A 3) 

C = 5; % Acoustic velocity of water (ft/msec) 

K1 = 22505; % Pmax 
A1 = 1.18; % Pmax 
K2 = 0.058; % Decay constant 
A2 = -0.185; % Decay constant 

% Plotting BC in 2D, TNT, 1000 lb @ 25ft: 
counter = 0; 
i = 1; 

for W = [1000]; % Equivalent charge weights (lb) 
for Dl= [25]; % Charge location depths (ft) 
counter = counter+1; 

A=zeros(50,700); 
for y = 1:51; 

for x = 1:701; 

R1 = sqrt((D1 - (y-l)) A 2 + (x-l) A 2); 

% Distance from charge to desired location (ft) 

R2 = sqrt((D1 + (y-l)) A 2 + (x-l) A 2); 

% Distance from image charge to desired location (ft) 
theta = K2 *(W A (1/3))*( ((W A (1/3))/Rl) A (A2)); 

% Decay Constant (msec) 

Pi = (Kl* (W A (1/3)/Rl) A (Al) ) * (exp (-(R2 -Rl)/(C*theta))); 

% Incident Pressure Wave (psi) 

Ph = Gamma*(y-1); % Hydrostatic Pressure at y (psi) 

Pr = (Kl*((W A (1/3)/R2) A (Al))); % Refelcted Pressure Wave (psi) 

F = Pi + Pa + Ph - Pr; % Upper Bulk Caviataion Boundary 

G1 = -Pi/ (C*theta)* (1+ ( ( (R2-2*D1*( (D1+ (y-1) )/R2) )/Rl)*(A2*R2/R1-A2- 

1) ) ) ; 

G2 = -(Al*Pi/Rl A 2)*(R2-2*D1*((D1+(y-1))/R2)); 

G3 = Gamma*((D1+(y-1))/R2) ; 

G4 = (A1/R2)*(Pi+Pa + Ph); 

G = G1 + G2 + G3 + G4; % Lower BC Boundary 
if F > 0.001; % Combine BC Boundaries 
if G < 0; 

_ A(y,x) = 1; _ 
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end 

end 

if G > 0; 

A(y,x) = 1; 

end 

end 

end 

temp(:,:,counter) = A; 
end 

charge=num2str(W); 
i=i+l; 

end 

figure4=figure( 'NumberTitle' , 'off', 'Name', 'Bulk Cavitation 

1000LB TNT Charge UNDEX at 25FT Depth', 'PaperOrientation' , 'Landscape') 

subplot1=subplot(2,1,1, 'Parent' ,figure4); 

box (subplot1, 'on' ); 

hold(subplotl, 'all' ); 

spy(temp(:,:,!)) 

title ('SPY View') 

xlabel (' Radius (ft)') 

ylabel (' Depth (ft) ') 

plot (figure4) 

%% 2-D Line 

gamma=0.03703; % seawater weight density (lb/in A 3) 

pa=14.7; % atmospheric pressure (psi) 

c=5; % acoustic velocity (ft/msec) 

charge_type = 2; 

D = 25; 

W = 1000; 

% TNT Data: 

chg_name=' TNT charge' ; 
depth=num2str(D); 
weight=num2str(W); 

Kl=22505; 

Al=l.18; 

K2=0.058; 

A2=-0.185; 
ub_data=[]; 
lb_data=[]; 

% Calculation of the upper boundary: 
for x=0:700 

for y=0:0.1:510 

rl=sqrt((D-y) A 2+x A 2); 
r2 = sqrt( (D+y) A 2+x A 2); 
theta=K2 *W A (1/3)*(W A (1/3)/rl) A A2 ; 

F1=K1*(W A (1/3)/rl) A Al*exp(-(r2-rl)/(c*theta)); 

F2=(gamma*y*12)-Kl*(W A (1/3)/r2) A Al+pa ; 

F=F1+F2 ; 
if F<=0 

ub_data=[ub_data;F x -y] ; 
break 

end 

end 

end 

% Calculation of the lower boundary: 
for x=0:700 


Region: 
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for y=0:0.1:510 

rl=sqrt((D-y) A 2+x A 2); 
r2 = sqrt( (D+y) A 2+x A 2); 
theta=K2 *W A (1/3)*(W A (1/3)/rl) A A2; 
pi=Kl*(W A (1/3)/rl) A Al*exp(-(r2-rl)/(c*theta)); 

Gl=—(pi/(c*theta))* (1+( ( (r2-(2*D*(D+y)/r2))/rl)*( ((A2*r2)/rl)-A2-1))); 
G2=-((Al*pi)/rl A 2)*(r2-2*D*((D+y)/r2)); 

G3=(gamma*12)*( (D+y)/r2); 

G4= (Al/r2)*(pi+pa+(gamma*y*12)); 

G=G1+G2+G3+G4; 
if G>=0 

lb_data=[lb_data; G x -y] ; 
break 

end 

end 

end 

% Plot upper and lower boundary data in 2D: 
subplot(2,1,2) 

plot(ub_data(:,2),ub_data(:,3), 'b' ,lb_data(:,2),lb_data(:,3), ' r' ) 

xlabel (' Radius [ft]') 

ylabel ( 'Depth [ft] ' ) 

title (['2D Line View']) 

grid on 

%% 3D Mesh & Peak Mesh: BC View of 1000 Ibf TNT Charge at 25ft 
figure 

r = ub_data(:,2) ; 

z = lb_data(:,3) ; 

w = ub_data(:,3) ; 
theta = 0:pi/2 0:2*pi; 

X = bsxfun(@times,r,cos(theta)); 

Y = bsxfun(@times,r,sin(theta)); 

Z = repmat(z,1,length(theta)) ; 

W = repmat(w,1,length(theta)); 
grid on 

% subplot(1,2,1) 
me s h(X, Y, Z, W) 
title(' 3D General View ' ) 
xlabel (' Radius (ft)') 
ylabel (' Radius (ft)') 
zlabel (' Depth [ft]') 

% subplot(1,2,2) 
figure 

meshc(X, Y, Z, W) 
title(' 3D Peak View') 

xlabel (' Radius (ft)') 
ylabel (' Radius (ft)' 
zlabel (' Depth [ft]') 
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APPENDIX E. SAMPLE GEMGRID CODES 


f or mat =0 0 0 0 3 

<X GRID LI N E S > 

D a t u ml n d e x = 1 
Dat um=0. 

N C e I I s = n a Si zeFi rstCel I =20. Si zeLastCel I =20. Rati o = n a Wi dt h =2 4 3 8 4 . #8 0 0 f t 

<END X GRID LI N E S > 

<Y GRID LI N E S > 

Dat u ml n d e x = 1 
Dat um=0. 

NC e I I s =1 Si zeFi rstCel I =na Si zeLastCel I =na Rati o =1. Wi d t h = 1. #2D Run 

# N C e I I s = n a SI zeFi rstCel I =20. Si zeLastCel I =20. Rat i o = n a Wi dt h=9144. #3 0 0 f t 

<END Y GRID LI N E S > 

<Z GRID LI N E S > 

Dat u ml n d e x = 1 
Dat um=- 1 6 0 0 2. 


NCel 1 s =n a 

Si zeFi 

rstCel 1 

=20. 

Si zeLastCel 1 

=20. 

Rati o = na 

Wi dt h =7 6 2. 

#2 5ft 

clay 

NCe 1 1 s = n a 
water 

Si zeFi 

rstCel 1 

=20. 

Si zeLastCel 1 

=20. 

Rati o = na 

Wi dth=15240. 

#5 0Of t 

NCel 1 s =na 
<END Z GRI D 

Si zeFi 
LI NES> 

rstCel 1 

=20. 

Si zeLastCel 1 

=20. 

Rati o = na 

Wi dt h =7 6 2. 

#2 5ft 

clay 


<OPTI ONS > 

#PI ot=Dys mas P 
<END OPTI ONS> 
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APPENDIX F. SAMPLE PREGEMINI CODES 


f or mat =0 0 0 1 7 
<0PTI ONS > 

Mode=st ar t # PreGemini mode (Start/Rezone/Overl ay/Reparti ti on) 

StartTi me=0. # Set initial time (REAL/precalc) 

Gravi ty=lg # Gravity, positive upward ( REAL/ 1 g) 

<END OPTI ONS> 

<GRI D> 

Coordi nat es =cyIi ndri cal 
x Ce I I s =1219 d x =0. x Da t u m=0 

y CeI I s =1 dy =0. y Dat u m=0 

z CeI 1 s =8 3 8 d z =0. z Da t u m=0 

Gri dFi I e = . / gri d/gri d. asc 

<END GRI D> 

<SUBGRI DS> 

XSubGri ds=4 XParti ti on=auto 
YSubGri ds=l YParti ti on=auto 
ZSubGri ds =3 ZPar ti ti on=auto 
<END SUBGRI DS> 

<BOUNDARY CONDI TI ONS > 

x mi n =wa I I xmax=free # wal I 2/wal I / freeng/free/REAL 

# y mi n =wa II ymax=f r ee 
z mi n =wa II zmax=free 
<END BOUNDARY CONDI TI ONS > 

<MATERI ALS> 

Ma t e r i a I I D = h e solid 
Ma t e r i a I I D = h e"ga s 
Ma t e r i a I I D=wa t e r 
Ma t e ri a I I D=a i r 
Ma t e r i a I I D=cI a y 
<END MATERI ALS> 

<BURN> 

Unbur ned=he solid Bur ned=he gas T i me =0. Ref P t = (0. , 0., - 1 5 2 4.) #5 0 f t z-deep 2D Run 
<END BURN> 

<HYDROSTATI C FI ELD> 

pRef =1. 0e+6 z R e f =0. zMax=max # Ref pressure, ref location, zMax (must be 

1st Ii ne) 

Materi al =ai r e i =e r e f z Mi n =0. # top state 

Materi al =water e i =e r e f zMi n=- 15240. # 2nd state 

Materi al =cl ay ei=eref zMi n=- 16002. # bottom state 

<END HYDROSTATI C FI ELD> 


< I NI TI A L S T AT E S > # state matl g frac rho e p aO u 

v w 

Statel D = he solid Materi al =he solid Rho=r hor ef e i =e r e f E OS v a r =( 0. ) # 

explosive state 

Statel D=wat er Materi a I =water Rho=r hor ef ei =e r ef EOSvar =( 0. ) # 

water state 

Statel D=ai r Materi a I =a i r Rho=r horef e i =e r e f EOSvar =( 0. ) # 

at mo sphere state 

Statel D=cl ay Materi al =cl ay Rho=r hor ef e i =e r e f EOSvar =( 0. ) 

<END I NITIAL STATES > 

<F L OWFI E L D> 

Opt i on=hydrost at i c 

Opti on=bal I st at e=he solid mass=11340. Ref P t =( 0. , 0., - 1 5 2 4.) #2 51 b @ 5 0ft 

<END F L OWF I ELD> 

<TEXT OUTPUT> 

i mi n = 1 i max = 1 j mi n = 1 j ma x = 1 kmi n = 1 kmax=l 
<END TEXT OUTPUT> 


Fi I e=hbx_lsol i d. mtl 
Fi I e = hbx_ 1, mtl 

Fi I e=ti I I water, mtl 
F i I e =ai r. mtl 
Fi I e=cl ay. mtl 


# Coordi nates: CARTESIAN, CYLINDRICAL, or SPHERICAL 

# File name for grid file specify (none/STRING) 


BE 
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APPENDIX G. SAMPLE GEMINI CODES 


f or mat =0 0 0 1 5 





<CAS E > 





Start Fi 1 e 

= . / pregemi ni / restart 000000 xxx.bin # Start file (use " xxx.bin 

ext ensi on 

for parallel 

runs) 




Ti 11 e =5 0 0 f t 

# 

Title (Limit 40 characters) 


Mode=f 1 ui 

d 


# Run mode (Fluid/Coupled) 


<END CAS E > 





<TERMl NATE> 





EndSt ep=9 9 9 9 9 9 


# Ter mi nation step (1 NT/none) 


EndTi me=. 

15 


# Ter mi nation time ( 1 NT/ none) 


d t Mi n =1. e 

- 12 


# Terminate if step size is less than thi 

value 

(sec) 





<T RAP > 





# x =6 0 8. 

44 i =1 z 

=- 1 2 3. 8 6 Var =p 

Di f =0. 01 


<END TRAP > 




<END TERMI NATE > 




<1 NTEGRATI 0N> 




C F L =. 4 5 



# CFL safety factor (0.9 for ID. 0.45 for 

2D and 3D) 

CFLI ni t =. 

05 


# Initial step CFL factor 


L i mi t e r =2 



# Li mi ter: 0.0- 2.0 


LagEqual i ze =s ou 


# Equalize after Lagrange step 


RemapEqual i ze =on 


# Equalize after remap step 


Pr ot ec t =1 



# Protect ( =0 off, >0 on) 


<END 1 NTEGRATI ON> 




<CELL HI STORY> 




i =1 

j =0. k =- 

1 5 2 4. 

Icharge center at 50ft z-deep 


x =3 0 4 8. 

y=0. 

z=- 25. 

#at x =2 5c m z-dee p 


x =4 5 7 2. 

y=0. 

z=- 25. 



x =6 0 9 6. 

y=0. 

z=- 25. 

#plus every 100ft to 80 0ft 


x =7 6 2 0. 

y=0. 

z=- 25. 



x =914 4. 

y=0. 

z=- 25. 

fexcept at 150, 200, 250, 

3 0 0, 3 5 0, 

4 0 0 f t 





x =1 0 6 6 8. 

y=o. 

z =- 25. 



x =1 2 19 2. 

y=0. 

z=- 25. 



x =1 5 2 4 0. 

y=o. 

z=- 25. 



x =1 8 2 8 8. 

y=o. 

z=- 25. 



x =2 1 3 3 6. 

y=o. 

z=- 25. 



x =2 4 3 8 4. 

y=0. 

z=- 25. 



x =3 0 4 8. 

y=o. 

z =- 1 5 2 4. 

#a t x =5 Of t z-deep 


x =4 5 7 2. 

y=o. 

z =- 1 5 2 4. 



x =6 0 9 6. 

y=o. 

z =- 1 5 2 4. 

#p 1 us every 100ft to 8 0 0ft 


x =7 6 2 0. 

y=o. 

z =- 1 5 2 4. 



x =914 4. 

y=o. 

z =- 1 5 2 4. 



x =1 0 6 6 8. 

y=o. 

z =- 1 5 2 4. 



x =1 2 19 2. 

y=o. 

z =- 1 5 2 4. 



x =1 5 2 4 0. 

y=o. 

z =- 1 5 2 4. 



x =1 8 2 8 8. 

y=o. 

z =- 1 5 2 4. 



x =2 1 3 3 6. 

y=o. 

z =- 1 5 2 4. 



x =2 4 3 8 4. 

y=0. 

z =- 1 5 2 4. 



x =3 0 4 8. 

y=o. 

z =- 213 4. 

#a t x =7 Of t z-deep 


x =4 5 7 2. 

y=o. 

z =- 2134. 



x =6 0 9 6. 

y=o. 

z =- 213 4. 

#pl us every 100ft to 8 0 0ft 


x =7 6 2 0. 

y=o. 

z =- 2134. 



x =914 4. 

y=o. 

z =- 213 4. 



x =1 0 6 6 8. 

y=o. 

z =- 2134. 



x =1 2 19 2. 

y=o. 

z =- 213 4. 



x =1 5 2 4 0. 

y=o. 

z =- 213 4. 



x =1 8 2 8 8. 

y=o. 

z =- 213 4. 



x =2 1 3 3 6. 

y=o. 

z=- 2134. 



x =2 4 3 8 4. 

y=0. 

z =- 213 4. 



<END CELL HI STORY> 




<CONTOUR P L OTS > 




St epFr eq = 

3 0 0 0 


# Plot step frequency (INT/none/table) 


T i me F r e q = 

0. 001 


# Plot time frequency (REAL/none/table) 
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<TARGET PLOT STEP5> 

3 0 0 0 

<END TARGET PLOT STEP S> 

<TARGET PLOT Tl MES> 

0 . 0 0 0 1 
0 . 001 
0 . 025 
0 . 030 
0 . 035 
0 . 040 
0.045 
0 . 050 
0 . 055 
0 . 060 
0 . 065 
0 . 070 
0 . 075 
0 . 080 
0 . 085 
0 . 090 
0 . 095 
0 . 100 
0 . 105 
0 . 110 
0 . 120 
0 . 130 
0.140 
0 . 150 

<END TARGET PLOT TI ME S > 

<END CONTOUR PLOTS> 


<RESTART > 

St epFreq=1000 
SaveFr eq=100 

use 2 to delete every other file) 
<END RESTART> 


# Restart file output step frequency (INT/none) 

# Retention interval (ex: use 1 to save every file, 


<TEXT OUTPUT> 

St epF r eq =500 
Ti meFr eq=1000. 

i mi n =1 i max = 1 j mi n = 1 j max = 1 kmin- 
<END TEXT OUTPUT> 


# Output step frequency (INT/none) 

# Output time frequency (INT/none) 
1 kmax=l 


# < C O U P L E D > 

#Un i t s =c m- g - s 
#Bac k Pr es s u r e =0. e +6 
#LoadFai I edSWI =on 
#Body Or i gi n =( 0. , 0. , 
I c m) 

#BodyXAx i s =( 1. , 0., 
Fluid coordinates (cm) 
#Bod y YAx i s =( 0. , 1., 
Fluid coordinates (cm) 
#Bod y ZAx i s =( 0. , 0. , 
Fluid coordinates (cm) 
#<END COUP L E D> 


# Back pressure inside SWI bodies (d/cm / '2) 

# Load failed SWI elements (on/off) 

# Location of Body coordinate system in Fluid mesh 


s=( 1. , 0. , 0. ) 
nates (cm) 

# 

Vector 

along 

Body 

coordi 

n a t e 

system x 

axis i 

s=( 0. , 1. , 0. ) 
nates (cm) 

# 

Vector 

along 

Body 

coordi 

n a t e 

system y 

axis i 

s=( 0. , 0. , 1. ) 
nates (cm) 

# 

Vector 

along 

Body 

coordi 

nat e 

system z 

axis i 


#<E L E ME NT HI STORY > 

#1 

#2 

#3 

#<END ELEMENT HI STORY > 


#<NODE HI STORY> 

#1 

#2 

#3 

#4 

#<END NODE HI STORY> 








APPENDIX H. SAMPLE GEMFIELD CODES 


For mat =0 0 0 0 4 
<Ca s e > 

Di r Out =. / # Subdirectory for output files. (Required. Make optional?) 

# Use for current working directory. 

Series =xyz # Series: Output file identifier (1-10 characters or 

( Opti onal ) 

F i I e F o r ma t =Dy s ma s P # DysmasP (latest format), DysmasP2010 

or DysmasP 2 0 0 7 
# 

# F i I e F o r ma t =T e c p I ot Fi I eSpi i t=Si ngl e # Tecplot (latest format), TecpiotlO, 

or Tec pi ot 9 

# F i I e F o r ma t =Tec p I ot Fi I eSpi i t = B y V a r # Tecplot (latest format), TecpiotlO, 

or Tec pi ot 9 

# Fi I eFormat=Tecpl ot Fi I eSpi i t=ByTi me # Tecplot (latest format), TecpiotlO, 

or Tec pI ot 9 

# Fi I eFormat=Tecpl ot Fi I eSpi i t=ByVarAndTi me # Tecplot (latest format), TecpiotlO, 

or Tec pi ot 9 

# Fi I e F o r mat =vt kxml Fi I eType=bi nary # v t k x ml (ascii or binary) or 

vt kl egacy (ascii only) 

# Fi I e F o r mat =vt kxml Fi I eType=asci i # v t k x ml (ascii or binary) or 

vtkl egacy (ascii only) 

# Fi I eFormat=vtkl egacy # vtkxml (ascii or binary) or 

vt kl egacy (ascii only) 

<End Case> 


#N o t e: all vars are optional 
<T I ME > 

tBeg=- 1. 000 t E n d =. 2 5 # t Beg, t End: Time window for output (sec) (before scaling 

or offset) 

# Default: Entire s i mu I a t i o n t i me 

tOffset=0. 00 tScal e=l. 00 # tOffset, tScale: Offset and scaling for time (Default: 

t Of f s et =0, t Sc a I e =1) 

# t < = = tScal e*( tOffset+t) 

MaxTi meRec or ds =5 0 0 0 0 # Max number of time records to generate. 

#N o t e: These target times are in ADDITION to any data at times within the time 
window (t Be g, t E n d) 

fspecified above. If you want only the target times then set t B e g and t E n d negative 
(or very large). 

<TI ME STEP TARGETS> 

#1815 
#6 0 0 0 

<END TI ME STEP TARGETS > 

<END TI ME> 


<OPTI ONS > 

S h o wF a i I ed =o n # Show failed body elements (on or off, default is on) 

# (Use "on" to preserve correct element numbering in Tecplot) 
Verbosity = l # Amount of screen output ( Def a u I t =1) 

<END OPTI ONS > 


# At least one plot variable is required 

# Default units are cgs 
<F L UI D VARI AB L E S > 

##T h e s e vars are only available as cell-averaged quantities 

V a r = p Uni t s = p s i #Units: psi, ksi, Pa, MPa, bar or scaling factor 

V a r = u Uni t s =1. #Units: ml s, ft/s, in/s or scaling factor 

V a r = v Uni t s =1. #Units: ml s, ft/s, in/s or scaling factor 

Var=w Uni ts=ft/s #Units: ml s, ft/s, in/s or scaling factor 

V a r =e k i n #Units: J, MJ or scaling factor 

##T h e s e vars are available as c el I - aver aged ( Mat N u m= 0) or per material ( Ma t N u m>0) 

V a r = r Ma t N u m=0 Uni t s = 1. 0 #Units: kg/ m^ 3 or scaling factor 

V a r = r Ma t N u m= 1 Uni t s = 1. 0 #Units: kg/m A 3 or scaling factor 

V a r = r Ma t N u m=2 Uni t s =1. 0 #Units: kg/ m^ 3 or scaling factor 

V a r =e Ma t N u m=0 Uni t s = 1. 0 #Units: J, MJ or scaling factor 

V a r =e Ma t N u m= 1 Uni t s = 1. 0 #Units: J, MJ or scaling factor 

V a r =e Ma t N u m=2 Uni t s =1. 0 #Units: J, MJ or scaling factor 

Va r =ei nt Mat Nu m=0 Un i t s =1. 0 _ #Un i t s: | , M| or scaling factor _ 
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Var =ei nt Mat Num=l 
Var =ei nt Mat Num=2 


Uni t s =1. 0 
Uni t s =1. 0 


#Uni ts: 
#Uni ts: 


or scaling factor 
or scaling factor 


##T h e s e v a r s are only available per material ( Ma t l\l u m>0) 
Var =f Mat Num=l #No units 

Var =f Mat Num=2 #No units 

<END FLUI D V A RI A B L E S > 


<BODY NODE VARI AB L ES > 

# V a r = u 

# V a r = v 

# Va r =w 

<END BODY NODE VARI ABL ES > 


Uni t s =1. 
Uni t s =1. 


<BODY ELEMENT VARI AB L E S > 

# V a r = p Uni 

# Var=pPos Uni 

# Var =pNeg Uni 

< E N D BODY ELEMENT VARI A B L E S > 


m/ s , f t / s , 

ml s , f t / s , 

ml s , f t / s , 


or seal 
or seal 
or seal 


n g factor 
n g factor 
n g factor 


Pa, MPa, bar or seal 

Pa, MPa, bar or seal 

Pa, MPa, bar or seal 


ng factor 
n g factor 
ng factor 


#0pt i ona! , 
<SUBDOMAI N> 
i Mi n =1 
y Mi n =0. 
k Mi n =1 
Of f s et =( 0. 


If left unspecified, the entire domain will be output 


<END SUBDOMAI N> 


i Max=9999999 
y Ma x =0. 
kMax=9999999 
0 . , 0 . ) 


i De I t a =1 
j De I t a =1 
k De I t a =1 


##################################################### 

# The following is a Key to the Plottable Variables 
##################################################### 

# Variable: all owed Mat Num (0. . . 


# 

all = 

all point variables 

NA 

# 

u = 

velocity in r, x-di r 

0 

# 

v = 

velocity in y, t h e t a -di r 

0 

# 

w = 

velocity in z-dir 

0 

# 

P = 

pressure 

0 

# 

r = 

dens i t y 

a n y 

# 

e = 

total specific energy 

a n y 

# 

e i n t = 

internal energy 

a n y 

# 

e k i n = 

kinetic energy 

0 

# 

t = 

t e mp e r a t u r e 

0 

# 

f = 

v o 1 u me fraction 

>0 

# 

i = 

i mpu 1 s e intensity (t i me in sec) 

0 

# 

eos #= 

EOS variable number # 

0 

feature, 

it is picked automatically 


n 

# 

s t a t = 

status of cell (active/blocked) 

0 

# 

ma t = 

ma t e r i a 1 id 

NA 


(for his only) 
(only 1 ma t e r i a I 


al I owed to have this 


(currently only available for P-alpha EOS) 
(for field only) 


########################################################## 

# The following are Scaling Factors for some common units 
########################################################## 

# Pressure 


1. 4 5 0 3 7 7 e- 005 

#ps i 

1. 450377e-008 

#ks i 

0 . 1 

#Pa 

1. e- 007 

#MP a 

1. e- 006 

#b a r 

Force 


1. e- 005 

#N 

1. e- 008 

#k N 

2. 2 4 8 0 8 9 e- 006 

#1 bf 


# Acceleration 


0. 01 

#m/ s A 2 

0. 0328083 

#ft/s A 2 

0. 3937008 

#i n/ s A 2 

Vel oc i t y 

0. 01 

#m/ s 

0. 0328083 

#f t / s 

0. 3937008 

#i n / s 

Di stance 

0. 01 

#m 

0. 0328083 

#f t 

0. 3937008 

# i n 

Vo 1 u me 

1. e- 006 

#m A 

















# 

3. 531467e-005 

#f t 3 

# 

0. 06102374 

#i n A 3 

# 

De n si t y 


# 

1 0 0 0 #kq/ m A 3 

# 

Energy 


# 

1. e - 0 0 7 #j 


# 

1. e- 013 #MJ 


# 

Mass 


# 

0.001 #kg 


# 

Impulse intensity 

# 

0. 1 

# P a * s 

# 

1. e- 007 

#MP a * s 

# 

1. e- 006 

#b a r * s 

# 

1. 4 5 0 3 7 7 e- 005 

#ps i * s 

# 

1. 4 5 0 3 7 7 e- 002 

#psi * ms 
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APPENDIX I. SAMPLE GEMHIS CODES 


For mat =0 0 0 0 4 
<Ca s e > 

Di r Out =. / # Subdirectory for output files. (Required. Make optional?) 

# Use for current working directory. 

Seri e s ="" # Series: Output file identifier (1-10 characters or 

( Opti onal ) 

Fi I eFormat =dysmasp 
or DysmasP 2 0 0 7 
# 

# F i I e F o r ma t =Tec p I ot 
for mat), TecpiotlO, or 

# F i I eFor ma t =Tec pI ot 
for mat), TecpiotlO, or 

# F i I eFor ma t =Tec pI ot 
or Tec pI ot 9 

# F i I eFor ma t =Tec pI ot 
or Tec pi ot 9 

# F i I eFor ma t =Tec pI ot 
or Tec pi ot 9 

# 

# FileFormat=plain 
<End Case> 


# DysmasP (latest format), DysmasP2010 


Fi I eSpi i tPt=Si ngl e 
Tec pi ot 9 

Fi I eSpi i tPt=Si ngl e 
Tec pi ot 9 

Fi I eSpi i t Pt = B y V a r 

Fi I eSpi i t Pt = B y P t 

Fi I eSpi i t Pt =ByVarAndPt 


Fi I eSpi i tGbl =Si ngl e # Tecplot (latest 

Fi I eSpi i t Gbl = B y V a r # Tecplot (latest 

# Tecplot (latest format), TecpiotlO, 

# Tecplot (latest format), TecpiotlO, 

# Tecplot (latest format), TecpiotlO, 


Fi I eSpi i tPt=ByVarAndPt # 


#N o t e: all v a r s are optional 
<Ti me> 

t Beg =- 1. 0 0 0 t End =. 15 
or offset) 

t Of f s et =0. 00 t Sc a I e =1. 00 
t Of f s et =0, t Sc a I e =1) 

Ma xTi me Rec o r d s =5 0 0 0 0 0 
(defaul t =5 0 0 0 0 0) 

RemoveOver I a p=on 

overlap (on or off, default is 

< E n d T i me > 


# t Beg, t End: Time window for output (sec) (before scaling 

# Default: Entire s i mu I a t i o n t i me 

# tOffset, tScale: Offset and scaling for time (Default: 

# t < = = t Sc a I e* (t Of f s et +t) 

# Set max num of time records (used for array allocation) 

# Flag for removing oldest data if there is a time 
on) 


# Default units are 
< P o i n t V a r i a b I e s > 
##Thes e va r s are onI 
Va r =p 

V a r = u 

V a r = v 

V a r =w 

Va r =e ki n 

V a r = i 


cgs 


y available as 
Uni t s =1. 

Uni t s =1. 

Uni t s =1. 

Uni t s =1. 


##T h e s e 

vars 

are avail 

able 

as cell 

V a r = r 


Ma t Nu m=0 

Un 

t s =1. 

0 

V a r = r 


Ma t Nu m=l 

Un 

ts=l. 

0 

V a r = r 


Ma t Nu m=2 

Un 

t s = 1 . 

0 

V a r =e 


Ma t Nu m=0 

Un 

t s =1 . 

0 

V a r =e 


Ma t Nu m=l 

Un 

ts=l. 

0 

V a r =e 


Ma t Nu m=2 

Un 

t s =1 . 

0 

#Va r =e 


Ma t Nu m=3 

Uni t s =1 

. 0 

V a r =ei 

nt 

Ma t Nu m=0 

Un 

t s =1 . 

0 

V a r =ei 

nt 

Ma t Nu m=l 

Un 

t s = 1 . 

0 

Var =ei 

nt 

Ma t Nu m=2 

Un 

t s =1. 

0 


cel I- averaged quantities 

#Units: psi, ksi, Pa, MPa, bar or scaling factor 
#Units: mis, ft/s, in/s or scaling factor 

#Units: mis, ft/s, in/s or scaling factor 

#Units: mis, ft/s, in/s or scaling factor 

#Un i t s: J , MJ or scaling factor 
#No units 

averaged ( Mat N u m= 0) or per material ( Ma t N u m>0) 


#Un i t s 

kg/ m A 3 

or 

seal 

i ng 

factor 

#Un i t s 

kg/ m A 3 

or 

seal 

i nc 

factor 

#Un i t s 

kg/ m A 3 

or 

seal 

i nc 

factor 

#Uni t s 

J - Mj 

or 

seal i 

ng 

factor 

#Un i t s 

J , Mj 

or 

seal i 

ng 

factor 

#Un i t s 

J , Mj 

or 

seal i 

ng 

factor 

#Un i t s: I , Ml 

or 

seal 

i n g 

factor 

#Un i t s 

J , Mj 

or 

seal i 

ng 

factor 

#Un i t s 

J , Mj 

or 

seal i 

ng 

factor 

#Un i t s 

J - Mj 

or 

seal i 

ng 

factor 


##T h e s e vars are only available 
Va r =f Ma t N u m= 1 

Va r =f Ma t N u m=2 

< E n d Point Vari abl es> 


per ma t e r i a I (Mat Nu m>0) 
#No units 
#No units 


#################################### 
#The remaining sections are optional 
#################################### 
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<Opti o n s > 

Char geScal e = 1. # Scaling for specie (Energylnt, Energy, Mass, Volume, Radius) 

(def aul t =1. ) 

Verbosi ty=l # IVerbose: Screen output (0 = n o, 1 =s o me, 2=1 o t s) 

pRefl mp=- 5. 00e7 # Impulse intensity ref pressure (dynes/ c m^2) (if unset then use 

value from 1st record in file) 

<End Opti ons > 


<GI o b a I Va ri a bI e s > 

#N o t e: Ma 11\1 u m=0 is for total of all materials 
Var = r adi us Mat Nu m=l Uni t s =1. 

Var=vol cav Mat Num=3 Uni t s = 1. 

#f o r body/ bl ock fixed mesh option 

# Var =BodyCG-X 

# Var = B o d y C G - Y 

# Var = BodyCG- Z 

# Var =BodyVei - X 

# Var=BodyVel - Y 

# Var =BodyVel - Z 

# Var =BodyF - X 

# Var=BodyF-Y 

# Var =BodyF- Z 

< E n d Global V a r i a b I e s > 


#Un i t s: m, ft, in or scaling factor 
#C a v i t at i on vol ume in water 


< B o d y Global Var i abl es > 

Var =f orce - x 

Var =f orce - y 
Var =f orce - z 
Var =f orce 
#Var =i mpuI se 

< E n d Body Global V a r i abl e s > 

< B o d y Element V a r i abl e s > 

Var =p 
# V a r = i 

< E n d Body Element V a r i abl e s > 

<Body Node Vari abl es > 

Var =f orce- x 
Var =f orce - y 
Var =f orce - z 
Var =f orce 
Var =u 

V a r = v 

V a r =w 
Var =v eI 
Var =x 

V a r =y 

V a r =z 

<End Body Node Vari abl es > 


##################################################### 

# The following is a Key to the Global Variables 
##################################################### 

# Variable: all owed Mat Num (0. . 

# a I I all global variables NA 

# dt = Gemini time step size NA 

# ma s s = ma s s any 

# energy = total energy any 

# radius = equivalent radius >0 

# vol = volume >0 

# volcav = cavitated volume >0 

#BodyVel-[xyz] NA (forbo 

# BodyCG- [xyz] NA (forbo 

# BodyF- [xyz] NA (forbo 


(for body/block f 
(for body/ bl oc k f 
(for body/ bl oc k f 


x e d me s h opt 
xed mes h opt 
xed mes h opt 


##################################################### 
# The following is a Key to the Point Variables 
##################################################### 


# V a r i a b 


a 1 

1 

poi nt 

var i abl es 

v e 1 

OC 

I t y I 

n r, x - d i r 

v e 1 

oc 

ity I 

n y, theta- 

vel 

oc 

Ity i 

n z - d i r 


al I owed Mat Num (0. . . NMat) 
NA 
0 
0 
0 









# 

p = 

pressure 

0 

# 

r = 

dens i t y 

a n y 

# 

e = 

total specific energy 

a n y 

# 

e i n t = 

internal energy 

a n y 

# 

e k i n = 

kinetic energy 

0 

# 

t = 

t e mp e r a t u r e 

0 

# 

f = 

volume fraction 

>0 


# eos # 
feature 

# 

# s t a t = 

# mat = 


i mp u I s e i i 
= EOS v a r i r 
it is pi i 

status of 
ma t e r i a I i 


n t e n s i t y (t i me in sec) 
able number # 
eked automatically 

cell (act i v e / bl ocked) 


(for his only) 

(only 1 material allowed to have this 

(currently only available for P-alpha EOS) 
(for field only) 


##################################################### 

# The following is a Key to the Body Global Variables 
##################################################### 

# force-x = Total body force in r,x-dir (sumof 

# force-y = Total body force in y,t h e t a -di r (sumof 

# force-z = Total body force in z-dir (sumof 

# force = Total body force magnitude (sum of 

# i mp u I s e = I mp u I s e (t i me in sec) 

##################################################### 

# The following is a Key to the Body Element Variables 
##################################################### 

# p = Element pressure (net) 

# i = I mp u I s e intensity (t i me in sec) 

# ppos = Element pressure on "positive" side (normal 
interface element 

# pneg = Element pressure on "negative" side (normal 
coupling interface element 

##################################################### 

# The following is a Key to the Body Node Variables 
##################################################### 


node forces) 
node forces) 
node forces) 
node forces) 


points toward " ey 
points away from 


of coupl i ng 


# force-x 

# force-y 

# force-z 

# force 

# u = N 

# v = N 

# w = N 

# v e I = N 

# x = N 

# y = N 

# z = N 


Node force 
Node force 
Node force 


n r, x - d i r 
n y, t het a-di r 
n z-dir 


Node force magnitude 


Node vel oc 
Node vel oc 
Node vel oc 
Node vel oc 
Node posit 
Node posit 


Node pos 


t y in r,x - d i r 
t y in y, t het a- di r 
t y in z-dir 
t y ma g n i t u d e 
on ( r, x- d i r ) 
on (y, t het a-di r) 
on (z-dir) 


########################################################## 

# The following are Scaling Factors for some common units 
########################################################## 

# Pressure 

# 1. 450377e- 005 #psi 

# 1. 450377e- 008 #ksi 

# 0. 1 #Pa 

# 1. e- 0 0 7 #MPa 

# l.e- 0 0 6 #b a r 

# Force 

# 1. e- 0 0 5 #N 

# l.e- 0 0 8 #k N 

# 2. 248089e- 006 #1 bf 

# Acceleration 

# 0.01 #m/s / '2 

# 0.0 3 2 8 0 8 3 #f t / s A 2 

# 0.3 9 3 7 0 0 8 #in/s / '2 

# Vel oci t y 

# 0.01 #m/s 

# 0.0 3 2 8 0 8 3 #ft/s 

# 0.3 9 3 7 0 0 8 #in/s 

# Distance 

# 0.01 #m 

# 0.0 3 2 8 0 8 3 #ft 

# 0.3 9 3 7 0 0 8 #i n 

# Vo I u me 

# l.e- 0 0 6 #m"3 

# 3. 5 3 1 4 6 7 e - 0 0 5 #f t "3 

# 0.0 6 1 0 2 3 7 4 #i n A 3 

# Density 

# 1 0 0 0 tkg/rrn 


Bj 












APPENDIX J. SAMPLE MATERIAL FILE 


f or mat =0 0 0 0 2 
<0PTI ONS > 

title="sand eos using mi egruni sen eos" 
eosType=mi eg 
<END OPTI ONS> 

<REFERENCE> 
r hoRef =2. 023 
ei Ref =0. 
c Ref =2 6 0 0 0 0. 0 
<END REFERENCE> 

<LI Ml TS> 

r hoMi n = 1. e- 02 
ei Mi n=- 9. 9e +99 
p Mi n =2. e +4 
<END LI Ml TS > 

<EOS VARS> 
g a mma 0 = . 9 7 
S =1.86 
r ho0 =2. 023 
ei 0=0. 0 

p 0 =0. 0 

c 0 2 =4. 0 4 8 e +10 
<END EOS VARS> 
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